Numerical Methods Course Notes
Version 0.11
(UCSD Math 174, Fall 2004)
Steven E. Pav 1
October 13, 2005
department of Mathematics, MC0112, University of California at San Diego, La Jolla, CA 92093-0112.
<spav@ucsd. edu> This document is Copyright © 2004 Steven E. Pav. Permission is granted to copy,
distribute and/or modify this document under the terms of the GNU Free Documentation License, Version
1.2 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front-
Cover Texts, and no Back-Cover Texts. A copy of the license is included in the section entitled "GNU Free
Documentation License" .
2
Preface
These notes were originally prepared during Fall quarter 2003 for UCSD Math 174, Numerical
Methods. In writing these notes, it was not my intention to add to the glut of Numerical Analysis
texts; they were designed to complement the course text, Numerical Mathematics and Computing,
Fourth edition, by Cheney and Kincaid [7]. As such, these notes follow the conventions of that
text fairly closely If you are at all serious about pursuing study of Numerical Analysis, you should
consider acquiring that text, or any one of a number of other fine texts by e.g., Epperson, Hamming,
etc. [3, 4, 5].
Figure 1: The chapter dependency of this text, though some dependencies are weak.
Special thanks go to the students of Math 174, 2003-2004, who suffered through early versions
of these notes, which were riddled with (more) errors.
Revision History
0.0 Transcription of course notes for Math 174, Fall 2003.
0.1 As used in Math 174, Fall 2004.
0.11 Added material on functional analysis and Orthogonal Least Squares.
Todo
• More homework questions and example problems.
• Chapter on optimization.
• Chapters on basic finite difference and finite element methods?
• Section on root finding for functions of more than one variable.
i
ii
Contents
Preface i
1 Introduction 1
1.1 Taylor's Theorem 1
1.2 Loss of Significance 4
1.3 Vector Spaces, Inner Products, Norms 5
1.3.1 Vector Space 5
1.3.2 Inner Products 6
1.3.3 Norms 7
1.4 Eigenvalues 8
1.4.1 Matrix Norms 8
Exercises 10
2 A "Crash" Course in octave/Matlab 13
2.1 Getting Started 13
2.2 Useful Commands 17
2.3 Programming and Control 19
2.3.1 Logical Forks and Control 20
2.4 Plotting 21
Exercises 24
3 Solving Linear Systems 25
3.1 Gaussian Elimination with Naive Pivoting 25
3.1.1 Elementary Row Operations 25
3.1.2 Algorithm Terminology 27
3.1.3 Algorithm Problems 28
3.2 Pivoting Strategies for Gaussian Elimination 29
3.2.1 Scaled Partial Pivoting 30
3.2.2 An Example 31
3.2.3 Another Example and A Real Algorithm 32
3.3 LU Factorization 33
3.3.1 An Example 34
3.3.2 Using LU Factorizations 36
3.3.3 Some Theory 37
3.3.4 Computing Inverses 37
3.4 Iterative Solutions 37
3.4.1 An Operation Count for Gaussian Elimination 37
iii
iv CONTENTS
3.4.2 Dividing by Multiplying 39
3.4.3 Impossible Iteration 40
3.4.4 Richardson Iteration 40
3.4.5 Jacobi Iteration 41
3.4.6 Gauss Seidel Iteration 42
3.4.7 Error Analysis 42
3.4.8 A Free Lunch? 44
Exercises 46
4 Finding Roots 49
4.1 Bisection 49
4.1.1 Modifications 49
4.1.2 Convergence 50
4.2 Newton's Method 50
4.2.1 Implementation 51
4.2.2 Problems 53
4.2.3 Convergence 53
4.2.4 Using Newton's Method 54
4.3 Secant Method 55
4.3.1 Problems 57
4.3.2 Convergence 58
Exercises 60
5 Interpolation 63
5.1 Polynomial Interpolation 63
5.1.1 Lagranges Method 63
5.1.2 Newton's Method 64
5.1.3 Newton's Nested Form 66
5.1.4 Divided Differences 67
5.2 Errors in Polynomial Interpolation 69
5.2.1 Interpolation Error Theorem 69
5.2.2 Interpolation Error for Equally Spaced Nodes 73
Exercises 75
6 Spline Interpolation 79
6.1 First and Second Degree Splines 79
6.1.1 First Degree Spline Accuracy 80
6.1.2 Second Degree Splines 81
6.1.3 Computing Second Degree Splines 82
6.2 (Natural) Cubic Splines 82
6.2.1 Why Natural Cubic Splines? 83
6.2.2 Computing Cubic Splines 84
6.3 B Splines 84
Exercises 87
CONTENTS v
7 Approximating Derivatives 89
7.1 Finite Differences 89
7.1.1 Approximating the Second Derivative 91
7.2 Richardson Extrapolation 91
7.2.1 Abstracting Richardson's Method 92
7.2.2 Using Richardson Extrapolation 93
Exercises 95
8 Integrals and Quadrature 97
8.1 The Definite Integral 97
8.1.1 Upper and Lower Sums 97
8.1.2 Approximating the Integral 99
8.1.3 Simple and Composite Rules 99
8.2 Trapezoidal Rule 100
8.2.1 How Good is the Composite Trapezoidal Rule? 101
8.2.2 Using the Error Bound 103
8.3 Romberg Algorithm 104
8.3.1 Recursive Trapezoidal Rule 106
8.4 Gaussian Quadrature 107
8.4.1 Determining Weights (Lagrange Polynomial Method) 107
8.4.2 Determining Weights (Method of Undetermined Coefficients) 108
8.4.3 Gaussian Nodes 109
8.4.4 Determining Gaussian Nodes 110
8.4.5 Reinventing the Wheel 112
Exercises 113
9 Least Squares 117
9.1 Least Squares 117
9.1.1 The Definition of Ordinary Least Squares 117
9.1.2 Linear Least Squares 118
9.1.3 Least Squares from Basis Functions 119
9.2 Orthonormal Bases 121
9.2.1 Alternatives to Normal Equations 122
9.2.2 Ordinary Least Squares in octave/Matlab 124
9.3 Orthogonal Least Squares 124
9.3.1 Computing the Orthogonal Least Squares Approximant 128
9.3.2 Principal Component Analysis 129
Exercises 132
10 Ordinary Differential Equations 135
10.1 Elementary Methods 135
10.1.1 Integration and 'Stepping' 136
10.1.2 Taylor's Series Methods 136
10.1.3 Euler's Method 137
10.1.4 Higher Order Methods 137
10.1.5 Errors 138
10.1.6 Stability 138
10.1.7 Backwards Euler's Method 141
vi CONTENTS
10.2 Runge-Kutta Methods 143
10.2.1 Taylor's Series Redux 144
10.2.2 Deriving the Runge-Kutta Methods 144
10.2.3 Examples 146
10.3 Systems of ODEs 146
10.3.1 Larger Systems 147
10.3.2 Recasting Single ODE Methods 147
10.3.3 It's Only Systems 149
10.3.4 It's Only Autonomous Systems 150
Exercises 154
A Old Exams 157
A.l First Midterm, Fall 2003 157
A. 2 Second Midterm, Fall 2003 158
A. 3 Final Exam, Fall 2003 159
A.4 First Midterm, Fall 2004 161
A. 5 Second Midterm, Fall 2004 162
A. 6 Final Exam, Fall 2004 164
B GNU Free Documentation License 167
1. APPLICABILITY AND DEFINITIONS 167
2. VERBATIM COPYING 168
3. COPYING IN QUANTITY 168
4. MODIFICATIONS 168
5. COMBINING DOCUMENTS 169
6. COLLECTIONS OF DOCUMENTS 169
7. AGGREGATION WITH INDEPENDENT WORKS 169
8. TRANSLATION 170
9. TERMINATION 170
10. FUTURE REVISIONS OF THIS LICENSE 170
ADDENDUM: How to use this License for your documents 170
Bibliography 171
Chapter 1
Introduction
1.1 Taylor's Theorem
Recall from calculus the Taylor's series for a function, f(x), expanded about some number, c, is
written as
f(x) ~ a + ai (x - c) + a 2 (x - c) 2 +
Here the symbol ~ is used to denote a "formal series," meaning that convergence is not guaranteed
in general. The constants a,{ are related to the function / and its derivatives evaluated at c. When
c = 0, this is a MacLaurin series.
For example we have the following Taylor's series (with c = 0):
2 S
sin (x) = x - — + — -.. .
x 2 x 4
cos(x) = l-- + ir -...
(1.1)
(1.2)
(1.3)
Theorem 1.1 (Taylor's Theorem). If /(x) has derivatives of order 0, 1, 2, . . . , n+1 on the closed
interval [a, 6] , then for any x and c in this interval
f(x) = f(k) (c) {x ~ c)k + /(n+1) ^ {x ~ c)n+1
k=0
k\
(n + 1)!
where ^ is some number between x and c, and f k (x) is the /c th derivative of / at x.
We will use this theorem again and again in this class. The main usage is to approximate
a function by the first few terms of its Taylor's series expansion; the theorem then tells us the
approximation is "as good" as the final term, also known as the error term. That is, we can make
the following manipulation:
1
2
CHAPTER 1. INTRODUCTION
k=0
/W (c) (x-c)<
fc=0
fc!
/W (
C ) {X-C) k (0 (X-C)
fc!
E
fc=0
/(n+D (g) (g _ c) n+l
(n + 1)!
1/^) (oiix- c r +i
(n + 1)!
n+1
(n+1)!
On the left hand side is the difference between f(x) and its approximation by Taylor's series.
We will then use our knowledge about /( n+1 ) (£) on the interval [a, b] to find some constant M such
that
/(*) - E
c x — c
fc=0
k\
|/(n+D (Q| |x-C
(n + 1)!
n+1
< M |x - c
n+1
Example Problem 1.2. Find an approximation for f(x) = sinx, expanded about c = 0, using
n = 3.
Solution: Solving for is fairly easy for this function. We find that
f(x) = sinx = sin(O) +
cos(O) x — sin(0) x 2 — cos(O) x 3 sin(£) x
1!
+
2!
+
3!
+
4!
so
x 3 sin(f ) x 4
sinx
sin(£) x 4
24
x 4
< — ,
- 24'
because |sin(£)| < 1.
Example Problem 1.3. Apply Taylor's Theorem for the case n = 1.
Solution: Taylor's Theorem for n = 1 states: Given a function, /(x) with a continuous derivative
on [a,b], then
m = f(c)+f(o(x-c)
for some £ between x, c when x, c are in [a, b].
This is the Mean Value Theorem. As a one-liner, the MVT says that at some time during a trip,
your velocity is the same as your average velocity for the trip. H
Example Problem 1.4. Apply Taylor's Theorem to expand f(x)
Solution: Simple calculus gives us
21x + 17 around c = 1.
\x
'(x
'(x
)(x
)(x
= x 3 - 21x 2 + 17,
= 3x 2 - 42x,
= 6x - 42,
= 6,
= 0.
1.1. TAYLOR'S THEOREM
3
with the last holding for k > 3. Evaluating these at c = 1 gives
f(x) = -3 + -39(x - 1) + y - '- +
6
Note there is no error term, since the higher order derivatives are identically zero. By carrying out
simple algebra, you will find that the above expansion is, in fact, the function f(x). H
There is an alternative form of Taylor's Theorem, in this case substituting x + h for x, and x
for c in the more general version. This gives
Theorem 1.5 (Taylor's Theorem, Alternative Form). If f(x) has derivatives of order 0, 1, ... , n+
1 on the closed interval [a,b], then for any x in this interval and any h such that x + h is in this
interval,
/(l + ft) ^Z!^ + /''7«)(f +1 ,
^ fc! (n+1)!
where £ is some number between x and x + h.
We generally apply this form of the theorem with h — > 0. This leads to a discussion on the
matter of Orders of Convergence. The following definition will suffice for this class
Definition 1.6. We say that a function f(h) is in the class O (h k ) (pronounced "big-Oh of h k ")
if there is some constant C such that
\f(h)\<C\h\ k
for all h "sufficiently small," i.e., smaller than some h* in absolute value.
For a function / G O (h k ) we sometimes write / = O (/i fc ) . We sometimes also write O (h h ),
meaning some function which is a member of this class.
Roughly speaking, through use of the "Big-O" function we can write an expression without
"sweating the small stuff." This can give us an intuitive understanding of how an approximation
works, without losing too many of the details.
Example 1.7. Consider the Taylor expansion of Inx:
w u . , (i/x) h (-i/* 2 ) /i 2 (2/e)
\n(x + h) = lnx+ v/ i 7 + ^ '-—I +
Letting x = 1, we have
ln(l + h) = h-^ + ^h\
Using the fact that £ is between 1 and 1 + h, as long as h is relatively small (say smaller than ^),
the term can be bounded by a constant, and thus
h 2
In (1 + ft) = h- — + 0(h 3 ).
Thus we say that h — 4^ is a O (/i 3 ) approximation to ln(l + h). For example
01 2
ln(l + 0.01) 0.009950331 w 0.00995 = 0.01 .
4
CHAPTER 1. INTRODUCTION
1.2 Loss of Significance
Generally speaking, a computer stores a number x as a mantissa and exponent, that is x = ±r x 10 fc ,
where r is a rational number of a given number of digits in [0.1, 1), and k is an integer in a certain
range.
The number of significant digits in r is usually determined by the user's input. Operations
on numbers stored in this way follow a "lowest common denominator" type of rule, i.e., precision
cannot be gained but can be lost. Thus for example if you add the two quantities 0.171717 and
0.51, then the result should only have two significant digits; the precision of the first measurement
is lost in the uncertainty of the second.
This is as it should be. However, a loss of significance can be incurred if two nearly equal
quantities are subtracted from one another. Thus if I were to direct my computer to subtract
0.177241 from 0.177589, the result would be .348 x 10~ 3 , and three significant digits have been lost.
This loss is called subtractive cancellation, and can often be avoided by rewriting the expression.
This will be made clearer by the examples below.
Errors can also occur when quantities of radically different magnitudes are summed. For exam-
ple 0.1234 + 5.6789 x 10 -20 might be rounded to 0.1234 by a system that keeps only 16 significant
digits. This may lead to unexpected results.
The usual strategies for rewriting subtractive expressions are completing the square, factoring,
or using the Taylor expansions, as the following examples illustrate.
Example Problem 1.8. Consider the stability of \Jx + 1 — 1 when x is near 0. Rewrite the
expression to rid it of subtractive cancellation.
Solution: Suppose that x = 1.2345678 x 10~ 5 . Then ^/xTT ~ 1.000006173. If your computer
(or calculator) can only keep 8 significant digits, this will be rounded to 1.0000062. When 1 is
subtracted, the result is 6.2 x 10 -6 . Thus 6 significant digits have been lost from the original.
To fix this, we rationalize the expression
r~ T a i — XT .v^TT+l x + 1-1
\Jx + 1 — 1 = yx + l — l-
'Vs+T+l Vx + 1 + 1 Vx + 1 + 1 '
This expression has no subtractions, and so is not subject to subtractive cancelling. When x =
1.2345678 x 10 -5 , this expression evaluates approximately as
1.2345678 x 10~ 5 fi
» 6.17281995 x 10~ 6
2.0000062
on a machine with 8 digits, there is no loss of precision. H
Note that nearly all modern computers and calculators store intermediate results of calculations
in higher precision formats. This minimizes, but does not eliminate, problems like those of the
previous example problem.
Example Problem 1.9. Write stable code to find the roots of the equation x 2 + bx + c = 0.
Solution: The usual quadratic formula is
-b ± yjb 2 - 4c
X ± = o
Supposing that & > c> 0, the expression in the square root might be rounded to b 2 , giving two
roots x + = 0, X- = —b. The latter root is nearly correct, while the former has no correct digits.
1.3. VECTOR SPACES, INNER PRODUCTS, NORMS
5
To correct this problem, multiply the numerator and denominator of x + by —b — \Jb 2 — 4c to get
2c
x +
-b - Vb 2 - 4c
Now if b » c > 0, this expression gives root x+ = —c/b, which is nearly correct. This leads to the
pair:
-b - Vb' 2 - 4c 2c
o > X +
X_ =
-6 - V6 2 - 4c
Note that the two roots are nearly reciprocals, and if x- is computed, x + can easily be computed
with little additional work. H
Example Problem 1.10. Rewrite e x — cosx to be stable when x is near 0.
Solution: Look at the Taylor's Series expansion for these functions:
cosx
tZ/^ tZJ^
1 + X+ ¥ + 3! + ¥ + 5f +
= x + x 2 + ||- + C (x 5 )
This expression has no subtractions, and so is not subject to subtractive cancelling. Note that we
propose calculating x + x 2 + x 3 /6 as an approximation of e x — cosx, which we cannot calculate
exactly anyway. Since we assume x is nearly zero, the approximate should be good. If x is very
close to zero, we may only have to take the first one or two terms. If x is not so close to zero, we
may need to take all three terms, or even more terms of the expansion; if x is far from zero we
should use some other technique. H
2 4 fi
X X X
~ 21 + IT ~ 6f +
1.3 Vector Spaces, Inner Products, Norms
We explore some of the basics of functional analysis which may be useful in this text.
1.3.1 Vector Space
A vector space is a collection of objects together with a binary operator which is defined over an
algebraic field. 1 The binary operator allows transparent algebraic manipulation of vectors.
Definition 1.11. A collection of vectors, V, with a binary addition operator, +, defined over V,
and a scalar multiply over the real field R, forms a vector space if
1. For each u, v G V, the sum u + v is a vector in V. (i.e., the space is "closed under addition.")
2. Addition is commutative: u + v = v + u for each u, v G V.
3. For each u G V, and each scalar a G R the scalar product au is a vector in V. (i.e., the space
is "closed under scalar multiplication.")
4. There is a zero vector G V such that for any u G V, Ou = 0, where is the zero of R.
5. For any u G V, lu = u, where 1 is the multiplicative identity of R.
6. For any u, v G V, and scalars a, (3 G R, both (a +k [3)u = au + [3u and a (u + v) = au + av
hold, where +m is addition in R. (i.e., scalar multiplication distributes in both ways.)
1 For the purposes of this text, this algebraic field will always be the real field, E, though in the real world, the
complex field, C, has some currency.
6
CHAPTER 1. INTRODUCTION
Example 1.12. The most familiar example of a vector space is R ra , which is the collection of
n-tuples of real numbers. That is u G R n is of the form [ui, u 2 , ■ ■ ■ , u n ] T , where iij £ 1 for
i = 1, 2, . . . , n. Addition and scalar multiplication over R are defined pairwise:
[ui,u 2 ,.. • ,n n ] T + [vi,v 2 , . . . ,v n ] T = [ui + v 1 ,u 2 + v 2 ,.. -,u n + v n ] T , and
a [ui,u 2 , . . . ,u n ] T = [aui,av,2, ■ ■ ■ ,au n ] T
Note that some authors distinguish between points in n-dimensional space and vectors in re-
dimensional space. We will use R n to refer to both of them, as in this text there is no need to
distinguish them symbolically.
Example 1.13. Let X C R fc be an closed, bounded set, and let H be the collection of all functions
from X to R. Then H forms a vector space under the "pointwise" defined addition and scalar
multiplication over R. That is, for u, v G H, u + v is the function in H defined by [u + v] (x) =
u(x) + v(x) for all x G X. And for u G H, a G R, era is the function in [era] (x) = a (u(x)).
Example 1.14. Let X C R fc be a closed, bounded set, and let Hq be the collection of all functions
from X to R that take the value zero on dX. Then H forms a vector space under the "pointwise"
defined addition and scalar multiplication over R. The only difference between proving Hq is a
vector space and the proof required for the previous example is in showing that Hq is indeed closed
under addition and scalar multiplication. This is simple because if x G dX, then [u + v] (x) =
u(x) + v(x) = + 0, and thus u + v has the property that it takes value on dX. Similarly for au.
This would not have worked if the functions of Hq were required to take some other value on dX,
like, say, 2 instead of 0.
Example 1.15. Let V n be the collection of all 'formal' polynomials of degree less than or equal to
n with coefficients from R. Then V n forms a vector space over R.
Example 1.16. The collection of all real-valued mxn matrices forms a vector space over the reals
with the usual scalar multiplication and matrix addition. This space is denoted as R mxn . Another
way of viewing this space is it is the space of linear functions which carry vectors of R n to vectors
of R m .
1.3.2 Inner Products
An inner product is a way of "multiplying" two vectors from a vector space together to get a scalar
from the same field the space is defined over (e.g., a real or a complex number). The inner product
should have the following properties:
Definition 1.17. For a vector space, V, defined over R, a binary function, (,), which takes two
vectors of V to R is an inner product if
1. It is symmetric: (v,u) = (u,v).
2. It is linear in both its arguments:
(au + [3v, w) = a (u, w) + [3 (v, w) and
(u, av + f3w) = a (u, v) + (3 (u, w) .
A binary function for which this holds is sometimes called a bilinear form.
3. It is positive: (v, v) > 0, with equality holding if and only if v is the zero vector of V.
1.3. VECTOR SPACES, INNER PRODUCTS, NORMS
7
Example 1.18. The most familiar example of an inner product is the L 2 (pronounced "L two")
inner product on the vector space R n . If u = [u±, 112, ■ ■ ■ , u n ] T , and v = [v 1, V2, ■ ■ ■ , v n ] T , then
letting
(u,v) 2 = y^ujVj
i
gives an inner product. This inner product is the usual vector calculus dot product and is sometimes
written as u ■ v or u v.
Example 1.19. Let H be the vector space of functions from X to R from Example 1.13. Then
for u,v G H, letting
{u,v) H = I u(x)v(x)dx,
Jx
gives an inner product. This inner product is like the "limit case" of the L 2 inner product on R n
as n goes to infinity.
1.3.3 Norms
A norm is a way of measuring the "length" of a vector:
Definition 1.20. A function ||-|| from a vector space, V, to R + is called a norm if
1. It obeys the triangle inequality: \\x + y\\ < \\x\\ + \\y\\ .
2. It scales positively: ||aa;|| = |a| , for scalar a.
3. It is positive: > 0, with equality only holding when x is the zero vector.
The easiest way of constructing a norm is on top of an inner product. If (, ) is an inner product
on vector space V, then letting
||ti|| = yj(u,u)
gives a norm on V. This is how we construct our most common norms:
Example 1.21. For vector x G R™, its L 2 norm is defined
n \ 2
This is constructed on top of the L 2 inner product.
Example 1.22. The L p norm on R n generalizes the L 2 norm, and is defined, for p > 0, as
/ n \ 1 /P
X
1= EN'
\i=l
Example 1.23. The L°° norm on R n is defined as
\ x \\™ = max .
00 i
The L°° norm is somehow the "limit" of the L p norm as p —>■ 00.
8
CHAPTER 1. INTRODUCTION
1.4 Eigenvalues
It is assumed the reader has some familiarity with linear algebra. We review the topic of eigenvalues.
Definition 1.24. A nonzero vector x is an eigenvector of a given matrix A, with corresponding
eigenvalue A if
Ax = Xx
Subtracting the right hand side from the left and gathering terms gives
(A — X\)x = 0.
Since x is assumed to be nonzero, the matrix A — Al must be singular. A matrix is singular if and
only if its determinant is zero. These steps are reversible, thus we claim A is an eigenvalue if and
only if
det (A - Al) = 0.
The left hand side can be expanded to a polynomial in A, of degree n where A is an n x n matrix. This
gives the so-called characteristic equation. Sometimes eigenvectors, -values are called characteristic
vectors,-values.
Example Problem 1.25. Find the eigenvalues of
" 1 1
4 -2
Solution: The eigenvalues are roots of
= det
1 - A
4
1
-2- A
This equation has roots Ai = — 3, A2 = 2.
(l-A)(-2-A)-4 = A 2 + A-6.
Example Problem 1.26. Find the eigenvalues of A 2 .
Solution: Let A be an eigenvalue of A, with corresponding eigenvector x. Then
A 2 x = A (Ax) = A (Ax) = AAa; = A 2 cc.
The eigenvalues of a matrix tell us, roughly, how the linear transform scales a given matrix; the
eigenvectors tell us which directions are "purely scaled." This will make more sense when we talk
about norms of vectors and matrices.
1.4.1 Matrix Norms
Given a norm ||-|| on the vector space, M. n , we can define the matrix norm "subordinate" to it, as
follows:
Definition 1.27. Given a norm ||-|| on M. n , we define the subordinate matrix norm on IR nxn by
IIAxll
A = max —r — TT-.
x^O \\X\\
1.4. EIGENVALUES
9
We will use the subordinate two-norm for matrices. From the definition of the subordinate
norm as a max, we conclude that if a; is a nonzero vector then
|Ax|
- < ||A|| 2 thus,
| J\tjC 1 1 2 1 1 1 1 2 1 1 *^ 1 1 2 *
Example 1.28. Strange but true: If A is the set of eigenvalues of A, then
Example Problem 1.29. Prove that
I AIL = max I A| .
uz AeA
|AB|| 2 < ||A|| 2 ||B|| 2 .
Solution:
• I I I A B X cy /\ IL II BX|L Mam i i i i
AB L =df max — — — < max - — ^- Tl — — = A L B L.
x^O X 2 x^O X L
10
CHAPTER 1. INTRODUCTION
Exercises
(1.1)
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
1.10)
1.11)
1.12)
1.13)
1.14)
1.15)
1.16)
1.17)
1.18)
1.19)
1.20)
1.21)
1.22)
1.23)
1.24)
Suppose / G O (h k ) . Show that / G O (h m ) for any m with < m < k. (Hint: Take
h* < 1.) Note this may appear counterintuitive, unless you remember that O (h k ) is a better
approximation than O (h m ) when m < k.
Suppose f £ O (h k ) , and g G O (ft™) . Show that fg £ O (h k+m ) .
Suppose / eO (/i fc ) , and g £ O (h m ) , with m < fc. Show that f + g £ O (h m ) .
Prove that /(/i) = -3/i 5 is in O (h 5 ).
Prove that f(h) = h 2 + 5/i 17 is in (/i 2 ).
Prove that f(h) = h 3 is not in (/i 4 ) (Hint: Proof by contradiction.)
Prove that sin(/i) is in (h).
Find a (h 3 ) approximation to sin/i.
Find a (/i 4 ) approximation to ln(l + /i). Compare the approximate value to the actual when
h = 0.1. How does this approximation compare to the O (/j 3 ) approximate from Example 1.7
for h = o.n
Suppose that / G O (h k ). Can you show that /' G O (h k ~ l )l
\/l to get rid of subtractive cancellation when x « 0.
yfx to get rid of subtractive cancellation when x is very large.
cosx of subtractive cancellation for x
5 ^ approximate.
cos 2 x of subtractive cancellation for x
Rewrite \fx + 1
Rewrite \Jx + 1
Use a Taylor's expansion to rid the expression 1
small. Use a O (x 5 ) approximate.
Use a Taylor's expansion to rid the expression 1 -
small. Use a (x 6 ) approximate.
Calculate cos(7r/2 + 0.001) to within 8 decimal places by using the Taylor's expansion.
Prove that if x is an eigenvector of A then ax is also an eigenvector of A, for the same
eigenvalue. Here a is a nonzero real number.
Prove, by induction, that if A is an eigenvalue of A then A fc is an eigenvalue of A fc for integer
k > 1. The base case was done in Example Problem 1.26.
Let B = Yli=o ai ^ l > wnere A = I. Prove that if A is an eigenvalue of A, then X^=o ai/ ^ ^ s
an eigenvalue of B. Thus for polynomial p(x), p(X) is an eigenvalue of p(A).
Suppose A is an invertible matrix with eigenvalue A. Prove that A -1 is an eigenvalue for
A" 1 .
Suppose that the eigenvalues of A are 1, 10, 100. Give the eigenvalues of B = 3A 3 — 4A 2 + I.
Show that B is singular.
Show that if || x \\ 2 = r, then x is on a sphere centered at the origin of radius r, in M. n .
If ||x|| 2 = 0, what does this say about vector xl
Letting x = [3 4 12] 1
What is the norm of
what is IIjcIL?
A =
1
1/2
1/3
1/n
(1.25) Show that ||A|| 2 = implies that A is the matrix of all zeros.
(1.26) Showthat || A" 1 1| 2 equals (1/| A
mini) > where \ m in is the smallest, in absolute value, eigenvalue
of A.
(1.27) Suppose there is some \i > such that, for a given A,
|At?|| 2 > fj,\\v\
2-
1.4. EIGENVALUES
11
for all vectors v.
(a) Show that fx < ||A|| 2 . (Should be very simple.)
(b) Show that A is nonsingular. (Recall: A is singular if there is some x / such that
Ax = 0.)
(c) Show that ||A _1 || 2 < (l//x) .
(1.28) If A is singular, is it necessarily the case that ||A|| 2 = 0?
(1.29) If ||A|| 2 > /i > does it follow that A is nonsingular?
(1.30) Towards proving the equality in Example Problem 1.28, prove that if A is the set of eigen-
values of A, then
IIAII > max |A| ,
AeA
where ||-|| is any subordinate matrix norm. The inequality in the other direction holds when
the norm is ||-|| 2 , but is difficult to prove.
CHAPTER 1. INTRODUCTION
Chapter 2
A "Crash" Course in octave/Matlab
2.1 Getting Started
Matlab is a software package that allows you to program the mathematics of an algorithm without
getting too bogged down in the details of data structures, pointers, and reinventing the wheel.
It also includes graphing capabilities, and there are numerous packages available for all kinds of
functions, enabling relatively high-level programming. Unfortunately, it also costs quite a bit of
money, which is why I recommend the free Matlab clone, octave, available under the GPL 1 , freely
downloadable from http : / /www . octave . org.
In a lame attempt to retain market share, Mathworks continues to tinker with Matlab to make
it noncompatible with octave; this has the side effect of obsoletizing old Matlab code. I will try
to focus on the intersection of the two systems, except where explicitly noted otherwise. What
follows, then, is an introduction to octave; Matlab users will have to make some changes.
You can find a number of octave/Matlab tutorials for free on the web; many of them are
certainly better than this one. A brief web search reveals the following excellent tutorials:
• http : //www . math .mtu . edu/~msgocken/intro/ intro . html
• http : //www. cyclismo . org/tutorial/matlab/vector .html
• http : //web . ew . usna . edu/~mecheng/DESIGN/CAD/MATLAB/usna . html
Matlab has some demo programs covering a number of topics- from the most basic functionality
to the more arcane toolboxes. In Matlab, simply type demo.
What follows is a lame demo for octave. Start up octave. You should get something like:
GNU Octave, version 2.1.44 (i686-pc-linux-gnu) .
Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003 John W. Eaton.
This is free software; see the source code for copying conditions.
There is ABSOLUTELY NO WARRANTY; not even for MERCH ANT I B I L I TY or
FITNESS FOR A PARTICULAR PURPOSE. For details, type 'warranty'.
Please contribute if you find this software useful.
For more information, visit http://www.octave.org/help-wanted.html
Report bugs to <bug-octave@bevo . che . wise . edu> .
octave : 1>
1 Gmi Public License. See http://www.gnu.org/copyleft/gpl.html.
13
14
CHAPTER 2. A "CRASH" COURSE IN OCT AVE /MATL AB
You now have a command line. The basic octavian data structure is a matrix; a scalar is a 1 x 1
matrix, a vector is an n x 1 matrix. Some simple matrix constructions are as follows:
octave :2> a = [12 3]
a =
12 3
octave :3> b = [5; 4; 3]
b =
5
4
3
octave :4> c = a'
c =
1
2
3
octave :5> d = 5*c - 2 * b
d =
-5
2
9
You should notice that octave "echoes" the lvalues it creates. This is either a feature or an
annoyance. It can be prevented by appending a semicolon at the end of a command. Thus the
previous becomes
octave :5> d = 5*c - 2 * b;
octave : 6>
For illustration purposes, I am leaving the semicolon off. To access an entry or entries of a
matrix, use parentheses. In the case where the variable is a vector, you only need give a single
index, as shown below; when the variable is a matrix, you need give both indices. You can also
give a range of indices, as in what follows.
WARNING: vectors and matrices in octave/Matlab are indexed starting from 1, and not
from 0, as is more common in modern programming languages. You are warned! Moreover, the
last index of a vector is denoted by the special symbol end.
octave :6> a(l) = 77
a =
77 2 3
octave :7> a(end) = -400
a =
77 2 -400
2.1. GETTING STARTED
15
octave :8> a(2:3) = [22 333]
a =
77 22 333
octave :9> M = diag(a)
M =
77
22
333
octave :10> M(2,l) = 14
M =
77
14 22
333
octave:ll> M(l:2,l:2) = [1 2;3 4]
M =
12
3 4
333
The command diag(v) returns a matrix with v as the diagonal, if v is a vector. diag(M) returns
the diagonal of matrix M as a vector.
The form c : d gives returns a row vector of the integers between a and d, as we will examine
later. First we look at matrix manipulation:
octave: 12> j = M * b
j =
13
31
999
octave : 13>
N = rand (3, 3)
N =
0.706307
0.911833
0.166880
0.027866 0.087402
0.624716 0.067067
0.769423 0.938714
octave : 14>
L =
L =
M + N
1.166880
3.706307
0.911833
2.027866 0.087402
4.624716 0.067067
0.769423 333.938714
octave :15> P = L * M
P =
16
CHAPTER 2. A "CRASH" COURSE IN OCT AVE /MATL AB
7.2505e+00
1.7580e+01
3.2201e+00
1.0445e+01
2.5911e+01
4.9014e+00
2.9105e+01
2.2333e+01
1.1120e+05
octave :16> P =
P =
1.1669e+00
1.1119e+01
0.0000e+00
L .* M
4.0557e+00
1.8499e+01
0.0000e+00
0.0000e+00
0.0000e+00
1.1120e+05
octave : 17> x = M \ b
x =
-6.0000000
5.5000000
0.0090090
octave :18> err = M * x - b
err =
Note the difference between L * M and L . * M; the former is matrix multiplication, the latter
is element by element multiplication, i.e.,
(L-*M). j = H j K id .
The command rand(m,n) gives an m x n matrix with each element "uniformly distributed" on
[0, 1]. For a zero mean normal distribution with unit variance, use randn(m,n) .
In line 17 we asked octave to solve the linear system
Mx = b,
by setting
x = M\b = M _1 b.
Note that you can construct matrices directly as you did vectors:
octave :19> B = [1 3 4 5; 2 -2 2 -2]
B =
13 4 5
2-2 2-2
You can also create row vectors as a sequence, either using the form c : d or the form c : e : d, which
give, respectively, c, c + 1, . . . , d, and c, c + e, . . . , d, (or something like it if e does not divide d — c)
as follows:
octave :20> z = 1:5
z =
1 2 3 4 5
2.2. USEFUL COMMANDS
17
octave:21> z = 5: (-1) : 1
z =
5 4 3 2 1
octave :22> z = 5:(-2):l
z =
5 3 1
octave :23> z = 2:3:11
z =
2 5 8 11
octave :24> z = 2:3:10
z =
2 5 8
Matrices and vectors can be constructed "blockwise." Blocks in the same row are separated by a
comma, those in the same column by a semicolon. Thus
octave :2> y=[2 7 9]
y =
2 7 9
octave :3> m = [z;y]
m =
2 5 8
2 7 9
octave :4> k = [(3:4) ' , m]
k =
3 2 5 8
4 2 7 9
2.2 Useful Commands
Here's a none too complete listing of useful commands in octave:
• help is the most useful command.
• floor (X) returns the largest integer not greater than X. If X is a vector or matrix, it computes
the floor element-wise. This behavior is common in octave: many functions which we normally
think of as applicable to scalars can be applied to matrices, with the result computed element-
wise.
• ceil(X) returns the smallest integer not less than X, computed element-wise.
• sin(X), cos(X), tan(X), atan(X), sqrt(X), returns the sine, cosine, tangent, arctangent,
square root of X, computed elementwise.
• exp(X) returns e x , elementwise.
• abs(X) returns |X| , elementwise.
CHAPTER 2. A "CRASH" COURSE IN OCT AVE /MATL AB
norm(X) returns the norm of X; if X is a vector, this is the L 2 norm:
/ \ 1/2
if X is a matrix, it is the matrix norm subordinate to the L 2 norm.
You can compute other norms with norm(X,p) where p is a number, to get the L p norm, or
with p one of Inf , -Inf , etc.
zeros (m,n) returns an m x n matrix of all zeros,
eye (m) returns the m x m identity matrix.
[m,n] = size (A) returns the number of rows, columns of A. Similarly the functions rows (A)
and columns (A) return the number of rows and columns, respectively.
length (v) returns the length of vector v, or the larger dimension if v is a matrix.
f ind(M) returns the indices of the nonzero elements of N. This may not seem helpful at first,
but it can be very useful for selecting subsets of data because the indices can be used for
selection. Thus, for example, in this code
octave :1> v = round(20*randn(400,3) ) ;
octave:2> selectv = v(f ind(v( : ,2) == 7),:)
we have selected the rows of v where the element in the second column equals 7. Now you
see why leading computer scientists refer to octave/Matlab as "semantically suspect." It is
a very useful language nonetheless, and you should try to learn its quirks rather than resist
them.
diag(v) returns the diagonal matrix with vector v as diagonal. diag(M) returns as a vector,
the diagonal of matrix v. Thus diag(diag(v) ) is v for vector v, but diag(diag(M) ) is the
diagonal part of matrix M.
toeplitz(v) returns the Toeplitz matrix associated with vector v. That is
" v(l) v(2) v(3)
v(2) v(l) v(2)
v(3) v(2) v(l)
toeplitz(v)
v(n - 1)
v(n - 2)
v(n) v(n-l) v(n-2) ••• v(l)
In the more general form, toeplitz (c,r) can be used to return a assymmetric Toeplitz
matrix.
A matrix which is banded on the cross diagonals is evidently called a "Hankel" matrix:
u(l) u(2) u(3) ••• u(n)
u(2) u(3) u(4) ••• v(2)
hankel(u,v)= u ( 3 ) u ( 4 ) u ( 5 ) v ( 3 )
u(n) v(2) v(3) ••• v(n)
eig(M) returns the eigenvalues of M. [V, LAMBDA] = eig(M) returns the eigenvectors, and
eigenvalues of M.
kron(M,N) returns the Kronecker product of the two matrices. This is a blcok construction
which returns a matrix where each block is an element of M as a scalar multiplied by the
whole matrix N.
flipud(N) flips the vector or matrix N so that its first row is last and vice versa. Similarly
fliplr(N) flips left/right.
2.3. PROGRAMMING AND CONTROL
19
2.3 Programming and Control
If you are going to do any serious programming in octave, you should keep your commands in a
file, octave loads commands from '.m' files. 2 If you have the following in a file called myfunc.m:
function [yl,y2] = myfunc(xl ,x2)
% comments start with a "/,'
% this function is useless, except as an example of functions.
% input :
% xl a number
% x2 another number
% output :
% yl some output
% y2 some output
yl = cos(xl) . * sin(x2);
y2 = norm(yl) ;
then you can call this function from octave, as follows:
octave :1> myfunc(2,3)
ans = -0.058727
octave :2> [a,b] = myfunc(2,3)
a = -0.058727
b = 0.058727
octave :3> [a,b] = myfunc([l 2 3 4] , [1 2 3 4] )
a =
0.45465 -0.37840 -0.13971 0.49468
b = 0.78366
Note this silly function will throw an error if xl and x2 are not of the same size.
It is recommended that you write your functions so that they can take scalar and vector input
where appropriate. For example, the octave builtin sine function can take a scalar and output a
scalar, or take a vector and output a vector which is, elementwise, the sine of the input. It is not
too difficult to write functions this way, it often only requires judicious use of . * multiplies instead
of * multiplies. For example, if the file myfunc.m were changed to read
yl = cos(xl) * sin(x2);
it could easily crash if xl and x2 were vectors of the same size because matrix multiplication is not
defined for an n x 1 matrix times another n x 1 matrix.
An .m file does not have to contain a function, it can merely contain some octave commands.
For example, putting the following into runner. m:
xl = rand (4, 3) ;
x2 = rand(size (xl) ) ;
[a,b] = myfunc(xl ,x2)
2 The 'm' stands for 'octave.'
20
CHAPTER 2. A "CRASH" COURSE IN OCT AVE /MATL AB
octave allows you to call this script without arguments:
octave :4> runner
a =
0.245936
0.246414
0.542728
0.257607
0.478054
0.186454
0.419457
0.378558
0.535323
0.206279
0.083917
0.768188
b = 1.3135
octave has to know where your .m file is. It will look in the directory from which it was called.
You can set this to something else with cd or chdir.
You can also use the octave builtin function f eval to evaluate a function by name. For example,
the following is a different way of calling myfunc.m:
octave :5> [a,b] = f eval ("myfunc" ,2,3)
a = -0.058727
b = 0.058727
In this form f eval seems like a way of using more keystrokes to get the same result. However, you
can pass a variable function name as well:
octave :6> fname = "myfunc"
fname = myfunc
octave :7> [a,b] = f eval (fname, 2, 3)
a = -0.058727
b = 0.058727
This allows you to effectively pass functions to other functions.
2.3.1 Logical Forks and Control
octave has the regular assortment of 'if-then-else' and 'for' and 'while' loops. These take the
following form:
if exprl
statements
elseif expr2
statements
elsif expr3
statements
else
statements
end
for var=vector
statements
end
while expr
2.4. PLOTTING
21
statements
end
Note that the word end is one of the most overloaded in octave/Matlab. It stands for the last
index of a vector of matrix, as well as the exit point for for loops, if statements, switches, etc. To
simplify debugging, it is also permissible to use endif to end an if statement, endfor to end a
for loop, etc..
The test expressions may use the logical conditionals: >, <, >=, <=, ==, ~=. Do not use the
assignment operator = as a conditional, as this can lead to disastrous results, as in C.
Here are some examples:
7oCompute the sign of x
if x >
s = 1;
elseif x ==
s = 0;
else
s = -1;
end
%a ' regular ' for loop
for i=l:10
sm = sm + i;
end
%an 'irregular' for loop
for i=[l 2 3 5 8 13 21 34]
f sm = f sm + i ;
end
while (sin(x) > 0)
x = x * pi;
end
2.4 Plotting
Plotting is one area in which there are some noticeable differences between octave and Matlab.
The commands and examples given herein are for octave, but the commands for Matlab are not
too different, octave ships its plotting commands to Gnuplot.
The main plot command is plot. You may also use semilogx, semilogy,loglog for 2D plots
with log axes, and contour and mesh for 3D plots. Use the help command to get the specific
syntax for each command. We present some examples:
n = 100;
X = pi .* ((l:n) ./ n);
Y = sin(X) ;
"/just plot Y
plot(Y) ;
22
CHAPTER 2. A "CRASH" COURSE IN OCT AVE /MATL AB
Zplot Y, but with the 'right' X axis labels
plot(X.Y) ;
W = sqrt(Y);
plot (W) ;
Zplot W, but with the 'right' X axis labels
plot(Y.W) ;
The output from these commands is seen in Figure 2.1. In particular, you should note the
difference between plotting a vector, as in Figure 2.1(c) versus plotting the same vector but with
the appropriate abscissa values, as in Figure 2.1(d).
(c) W = W (d) W versus Y
Figure 2.1: Four plots from octave.
Some magic commands are required to plot to a file. For octave, I recommend the following
magic formula, which replots the figure to a file:
%call the plot commands before this line
gset term postscript color;
2.4. PLOTTING
23
gset output "filename.ps";
replot ;
gset term xll;
gset output "/dev/null" ;
In Matlab, the commands are something like this:
%call the plot commands before this line
print (gcf , ' -deps ' , ' f ilename . eps ' ) ;
24
CHAPTER 2. A "CRASH" COURSE IN OCT AVE /MATL AB
Exercises
(2.1) What do the following pieces of octave/Matlab code accomplish?
(a) x = (0:40) ./ 40;
(b) a = 2;
b = 5;
x = a + (b-a) .* (0:40) ./ 40;
(c) x = a + (b-a) .* (0:40) ./ 40;
y = sin(x) ;
plot(x,y) ;
(2.2) Implement the naive quadratic formula to find the roots of x 2 + bx + c = 0, for real b, c. Your
code should return
-b ± Vfr 2 - 4c
2 '
Your m-file should have header line like:
function [xl,x2] = naivequad(b , c)
Test your code for (b, c) = (l x 10 15 , l) . Do you get a spurious root?
(2.3) Implement a robust quadratic formula (cf. Example Problem 1.9) to find the roots of x 2 +
bx + c = 0. Your m-file should have header line like:
function [xl,x2] = robustquad(b,c)
Test your code for (b, c) = (l x 10 15 , l) . Do you get a spurious root?
(2.4) Write octave/Matlab code to find a fixed point for the cosine, i.e., some x such that x =
cos(x). Do this as follows: pick some initial value xq, then let Xj+i = cos(xj) for i =
0,1,..., n. Pick n to be reasonably large, or choose some convergence criterion (i.e., ter-
minate if |xj+i — Xi\ < 1 x 1CP 10 ). Does your code always converge?
(2.5) Write code to implement the factorial function for integers:
function [nfact] = factorial (n)
where n factorial is equal to 1 ■ 2 ■ 3 ■ ■ • (n — 1) ■ n. Either use a 'for' loop, or write the function
to recursively call itself.
Chapter 3
Solving Linear Systems
A number of problems in numerical analysis can be reduced to, or approximated by, a system of
linear equations.
3.1 Gaussian Elimination with Naive Pivoting
Our goal is the automatic solution of systems of linear equations:
a n xi + 012^2 + 013X3 +
a 2 \X\ + 022^2 + 023^3 +
a 3i xi + 03222 + 033^3 +
a n \x\ + a n2 X2 + a n3 x 3 +
- (l\ n X n
+ a2nX n
3n X n
6l
b 2
b 3
&nnXn — b r ,
In these equations, the and bi are given real numbers. We also write this as
Ax = b,
where A is a matrix, whose element in the i th row and j th column is a^-, and b is a column vector,
whose z th entry is b^.
This gives the easier way of writing this equation:
Oil
ai2
a 13 ■ ■
Gin
Xi
' h '
021
022
023 • •
G2n
X 2
b 2
031
032
033 • •
• G3n
x 3
b 3
«n2
«n3 • •
"'Till
Xn
b n
3.1.1 Elementary Row Operations
You may remember that one way to solve linear equations is by applying elementary row operations
to a given equation of the system. For example, if we are trying to solve the given system of
equations, they should have the same solution as the following system:
25
26
CHAPTER 3.
SOLVING LINEAR SYSTEMS
fit 1
f7 1 o
"'In
A l
an
022
023 • •
02n
^2
031
033 • •
03n
^3
nan
na i2
na i3 ■ ■
K(li n
ubi
dnl
a n 2
«n3
2-n
where k is some given number which is not zero. It suffices to solve this system of linear
equations, as it has the same solution(s) as our original system. Multiplying a row of the system
by a nonzero constant is one of the elementary row operations.
The second elementary row operation is to replace a row by the sum of that row and a constant
times another. Thus, for example, the following system of equations has the same solution as the
original system:
a u
au
ai 3
ain
Xl
' h
an
022
«23
a2n
X2
032
«33
03n
%3
bz
a (i-l)2
«(i-l)3
a (i-l)n
X(i-l)
b(i-i)
an + fiaji
ai2 + (3aj 2
a i3 + (3a j% ■ ■
ain Pajn
x%
h + f3b j
a n i
a n 2
a n 3
ann
bn
We have replaced the i row by the i row plus (3 times the j row.
The third elementary row operation is to switch rows:
an
ai2
an • •
ain
Xl
' h '
au
032
a33 • •
a-3n
X2
b 3
021
a22
023 • •
a2n
X3
b 2
a n i
a n 2
a n 3 ■ ■
ann
Xn
b n
We have here switched the second and third rows. The purpose of this e.r.o. is mainly to make
things look nice.
Note that none of the e.r.o. 's change the structure of the solution vector x. For this reason, it is
customary to drop the solution vector entirely and to write the matrix A and the vector b together
in augmented form:
(
an
ai2
ai3 ■ ■
a In
bi
a2i
022
023 • •
«2n
b 2
a3i
032
<333 • •
03n
b 3
V
a n i
a n 2
«n3 - -
ann
b n
The idea of Gaussian Elimination is to use the elementary row operations to put a system into
upper triangular form then use back substitution. We'll give an example here:
3. 1 . GAUSSIAN ELIMINATION WITH NAIVE PIVOTING
27
Example Problem 3.1. Solve the set of linear equations:
x\ + x 2 - x 3 = 2
— 3xi — 4x 2 + 4x3 = — 7
2xi + lx 2 + lx 3 = 7
Solution: We start by rewriting in the augmented form:
1
-1
2
-3
-4
4
-7
2
1
1
7
We add 3 times the first row to the second, and —2 times the first row to the third to get:
-1
1
3
We now add —1 times the second row to the third row to get:
-1
1
2
The matrix is now in upper triangular form: there are no nonzero entries below the diagonal. This
corresponds to the set of equations:
x\ + x 2 - x 3 = 2
-X2 + x 3 = -1
2x 3 = 4
We now solve this by back substitution. Because the matrix is in upper triangular form, we can
solve X3 by looking only at the last equation; namely X3 = 2. However, once X3 is known, the second
equation involves only one unknown, x 2 , and can be solved only by x 2 = 3. Then the first equation
has only one unknown, and is solved by x\ = 1. H
All sorts of funny things can happen when you attempt Gaussian Elimination: it may turn
out that your system has no solution, or has a single solution (as above), or an infinite number
of solutions. We should expect that an algorithm for automatic solution of systems of equations
should detect these problems.
3.1.2 Algorithm Terminology
The method outlined above is fine for solving small systems. We should like to devise an algorithm
for doing the same thing which can be applied to large systems of equations. The algorithm will
take the system (in augmented form):
/
an
ai2
ai 3 • •
bi
\
a.21
022
a23 • •
0-2n
b 2
031
032
033 • •
03n
bz
V
O-nl
a n 2
a«3 • •
b n
J
28
CHAPTER 3. SOLVING LINEAR SYSTEMS
The algorithm then selects the first row as the pivot equation or pivot row, and the first element of
the first row, an is the pivot element. The algorithm then pivots on the pivot element to get the
system:
V o
an ai3
a
22
J 32
h n2
a 23
a
33
l fl3
a\ n
bi
\
a 2n
b' 3
Where
b'
b'n )
(2<i<n,l<j<n)
Effectively we are carrying out the e.r.o. of replacing the i th row by the i th row minus (f^-) ti
times
the first row. The quantity is the multiplier for the i th row.
Hereafter the algorithm will not alter the first row or first column of the system. Thus, the
algorithm could be written recursively. By pivoting on the second row, the algorithm then generates
the system:
/ an
an
Ol3 •
a\ n
bi
\
a 22
a 23 '
a 2n
b'2
a 33 '
b'l
V o o
l n3
In this case
% =
bi =
l 2j
b'
K )
(3 < i < n, 1 < j < n)
3.1.3 Algorithm Problems
The pivoting strategy we examined in this section is called 'naive' because a real algorithm is a bit
more complicated. The algorithm we have outlined is far too rigid-it always chooses to pivot on
the k th row during the k th step. This would be bad if the pivot element were zero; in this case all
the multipliers — are not defined.
Bad things can happen if a^k is merely small instead of zero. Consider the following example:
Example 3.2. Solve the system of equations given by the augmented form:
-0.0590 0.2372
0.1080 -0.4348
-0.3528
0.6452
Note that the exact solution of this system is x\ = 10, X2 = 1. Suppose, however, that the algorithm
uses only 4 significant figures for its calculations. The algorithm, naively, pivots on the first
equation. The multiplier for the second row is
ai08 ° -1.830508...,
-0.0590
which will be rounded to —1.831 by the algorithm.
3.2. PIVOTING STRATEGIES FOR GAUSSIAN ELIMINATION
29
The second entry in the matrix is replaced by
-0.4348 - (-1.831)(0.2372) = -0.4348 + 0.4343 = -0.0005,
where the arithmetic is rounded to four significant figures each time. There is some serious sub-
tractive cancellation going on here. We have lost three figures with this subtraction. The errors
get worse from here. Similarly, the second vector entry becomes:
0.6452 - (-1.831)(-0.3528) = 0.6452 - 0.6460 = -0.0008,
where, again, intermediate steps are rounded to four significant figures, and again there is subtrac-
tive cancelling. This puts the system in the form
-0.0590 0.2372
-0.0005
-0.3528
-0.0008
When the algorithm attempts back substitution, it gets the value
-0.0008
X2 = ^00005 =
This is a bit off from the actual value of 1. The algorithm now finds
xi = (-0.3528 - 0.2372 • 1.6) /-0.059 = (-0.3528 - 0.3795) /-0.059 = (-0.7323) /-0.059 = 12.41,
where each step has rounding to four significant figures. This is also a bit off.
3.2 Pivoting Strategies for Gaussian Elimination
Gaussian Elimination can fail when performed in the wrong order. If the algorithm selects a zero
pivot, the multipliers are undefined, which is no good. We also saw that a pivot small in magnitude
can cause failure. As here:
exi +x 2 = l
xi + x 2 = 2
The naive algorithm solves this as
l-X 2 1
Xi = =
e 1 — e
If e is very small, then ^ is enormous compared to both 1 and 2. With poor rounding, the algorithm
solves X2 as 1. Then it solves xi = 0. This is nearly correct for X2, but is an awful approximation
for xi. Note that this choice of xi,X2 satisfies the first equation, but not the second.
Now suppose the algorithm changed the order of the equations, then solved:
Xi + X2 = 2
exi + x 2 = 1
30
CHAPTER 3. SOLVING LINEAR SYSTEMS
The algorithm solves this as
1 - 2e
X2 = —
Xl = 2 - x 2
There's no problem with rounding here.
The problem is not the small entry per se: Suppose we use an e.r.o. to scale the first equation,
then use naive G.E.:
1 1
Xl + -Xl =
i €
xi + x 2 = 2
This is still solved as
Xl
and rounding is still a problem.
3.2.1 Scaled Partial Pivoting
The naive G.E. algorithm uses the rows 1,2,..., n-1 in order as pivot equations. As shown above,
this can cause errors. Better is to pivot first on row £\, then row £2, etc, until finally pivoting on
row £ n -i, for some permutation {£i}™ =1 of the integers 1, 2, . . . , n. The strategy of scaled partial
pivoting is to compute this permutation so that G.E. works well.
In light of our example, we want to pivot on an element which is not small compared to other
elements in its row. So our algorithm first determines "smallness" by calculating a scale, row-wise:
Si = max lajJ .
i<i<«
The scales are only computed once.
Then the first pivot, £1, is chosen to be the i such that
Si
is maximized. The algorithm pivots on row £1, producing a bunch of zeros in the first column. Note
that the algorithm should not rearrange the matrix-this takes too much work.
The second pivot, £2, is chosen to be the i such that
Si
is maximized, but without choosing £2 = l\. The algorithm pivots on row £2, producing a bunch of
zeros in the second column.
In the k th step £k is chosen to be the i not among £i,£2,..., £k~i such that
Si
2-±
e_
1 — X2
3.2. PIVOTING STRATEGIES FOR GAUSSIAN ELIMINATION
31
is maximized. The algorithm pivots on row £^, producing a bunch of zeros in the k column.
The slick way to implement this is to first set £i = i for i = 1,2, ... ,n. Then rearrange this
vector in a kind of "bubble sort": when you find the index that should be l\, swap them, i.e., find
the j such that £j should be the first pivot and switch the values of £\,£j.
Then at the fc th step, search only those indices in the tail of this vector: i.e., only among £j for
k < j < n, and perform a swap.
3.2.2 An Example
We present an example of using scaled partial pivoting with G.E. It's hard to come up with an
example where the numbers do not come out as ugly fractions. We'll look at a homework question.
1 2
-1
3
7
15 \
4
4
7
11
2
1
1
3
7
V 6
5
4
17
31 j
The scales are as follows: s\ = 7, s 2 = 7, S3 = 3, S4 = 17.
We pick i\. It should be the index which maximizes |aji| /sj. These values are:
2 4 2 6
7' 7'3' 17'
We pick £i=3, and pivot:
-2
2
4
8 ^
2
-2
1
-3
2
1
1
3
7
V
2
1
8
10 J
We pick £2- It should not be 3, and should be the index which maximizes (a^l /si- These values
are:
2 2 2
7' 7' 17'
We have a tie. In this case we pick the second row, i.e., £2 = 2. We pivot:
(
5
5
\
2
-2
1
-3
2
1
1
3
7
V
3
7
13
The matrix is in permuted upper triangular form. We could proceed, but would get a zero
multiplier, and no changes would occur.
If we did proceed we would have £3 = 4. Then l\ = 1. Our row permutation is 3, 2,4, 1. When
we do back substitution, we work in this order reversed on the rows, solving x±, then X3,X2,x±.
We get X4 = 1, so
^3 = ^ (13 -7*1) = 2
x 2 = - (-3- 1 * 1 + 2*2) =
xi = - (7 - 3*1-1*2-1*0) = !
32
CHAPTER 3. SOLVING LINEAR SYSTEMS
3.2.3 Another Example and A Real Algorithm
Sometimes we want to solve
Ax = b
for a number of different vectors b. It turns out we can run G.E. on the matrix A alone and come up
with all the multipliers, which can then be used multiple times on different vectors b. We illustrate
with an example:
/ 1 2 4 1 \ rr
4 2 12 2
2 12 3' 3
\ 1 3 2 1 / [4
The scale vector is s = [ 4 4 3 3] T .
Our scale choices are |, |, |, |. We choose £\ = 2, and swap £±,£2- In the places where there
would be zeros in the real matrix, we will put the multipliers. We will illustrate them here boxed:
Mr
Mi =
(
1
3
15
1 \
4
2
4
2
' 2
4
2
1
2
1
1
2
3
2
2
, i =
3
4
1
5
7
1
K
4
2
4
5 /
Our scale choices are §,§,§• We choose £2 = 4, and so swap £2, £4:
M 2 =
1
3
27
4
5
10
4
2
1
1
3
2
2
1
5
7
4
2
4
1 \
5
2
S 2
Our scale choices are |g, \ . We choose ^3 = 1, and so swap ^3,^4:
M,
V
27
10
1 \
5
2
11
9
1
2
2
4
1
3
Now suppose we had to solve the linear system for b = [ — 1 8 2 l] T -
We scale 6 by the multipliers in order: £\ = 2, so, we sweep through the first column of M3,
picking off the boxed numbers (your computer doesn't really have boxed variables), and scaling b
3.3. LU FACTORIZATION
33
appropriately:
This continues:
" -1 "
' -3 "
8
2
8
-2
1
_ -1 _
" -3 "
8
r 12 "
12 -i
5
8
5
8
-2
=>
-2
2
3
_ -1 _
-1
-1
27
1
12
\
10
5
5
1
2
8
7
4
17
2
?
2
3
-1
/
We then perform a permuted backwards substitution on the augmented system
/
4 2
This proceeds as
X4
£3
^2
Xi
-2 9 -6
~Tl7 ~ T7
10/12 1 -6
27 V T ~ 5 T7
2 / 1-6 7
I ' 8
-6
2— - x 3 - 2x 2
Fill in your own values here.
3.3 LU Factorization
We examined G.E. to solve the system
where A is a matrix:
Ax = b,
an «12 «13
«21 «22 0,23
«31 «32 ^33
a In
02n
«3n
Onl «n2 On3 '
We want to show that G.E. actually factors A into lower and upper triangular parts, that is A = LU,
where
1
•
••
21
1
•
••
31
^32
1 •
••
f"aV ^n2 ^n3 ' 1
We call this a LU Factorization of A.
U
U U U12 Ui 3
u 2 2 U23
u 33
Uln
U2n
U3n
34
CHAPTER 3. SOLVING LINEAR SYSTEMS
3.3.1 An Example
We consider solution of the following augmented form:
/ 2
1
1
3
7 \
4
4
7
11
6
5
4
17
31
V 2
-1
7
15 /
(3.2)
The naive G.E. reduces this to
1 2
1
1
3
7
\
2
-2
1
-3
3
7
13
V o
12
18
/
We are going to run the naive G.E., and see how it is a LU Factorization. Since this is the naive
version, we first pivot on the first row. Our multipliers are 2, 3, 1. We pivot to get
( 2
i
1
3
7
2
-2
1
-3
2
1
8
10
V o
-2
-1
4
8
/
Careful inspection shows that we've merely multiplied A and b by a lower triangular matrix Mi:
1 "
-1 1 _
The entries in the first column are the negative e.r.o. multipliers for each row. Thus after the first
pivot, it is like we are solving the system
MiAz; = Mib.
We pivot on the second row to get:
( 2
1
1
3
7
\
2
-2
1
-3
3
7
13
V o
-3
5
5
/
The multipliers are 1,-1. We can view this pivot as a multiplication by M2, with
M 2 =
1
1
-1
1
1
1
We are now solving
M 2 MiAx = M 2 Mib.
3.3. LU FACTORIZATION
35
We pivot on the third row, with a multiplier of —1. Thus we get
We have multiplied by M3 :
We are now solving
But we have an upper triangular form, that is, if we let
/ 2
1
_L
1
_L
O
7
\
2
-2
1
-3
3
7
13
V °
12
18 j
" 1
"
M 3
1
1
_
1
1 _
yi 3 M 2
Mi
Ace =
= M 3
M 2 M
lb
2
1
-2
3
3
1
7
12
Then we have
M 3 M 2 MiA= U,
A= (M 3 M 2 Mi) _1 U,
A
A
Mi"
LU.
• 1 M 2 " 1 M 3 - 1 U,
We are hoping that L is indeed lower triangular, and has ones on the diagonal. It turns out that
the inverse of each Mj matrix has a nice form (See Exercise (3.6)). We write them here:
1
"
" 1
"
" 1
2
1
1
1
3
1
1
1
1
1
1
-1
1
-1
1
1
2
3
1
1
-1
This is really crazy: the matrix L looks to be composed of ones on the diagonal and multipliers
under the diagonal.
Now we check to see if we made any mistakes:
LU
1
2
1
1
3
2
1
2
-2
1
3
1
1
3
7
1
-1
-1 1 _
_
12
2
1
1 3 "
4
6
4
5
7
4 17
= A.
2
-1
7
36
CHAPTER 3. SOLVING LINEAR SYSTEMS
3.3.2 Using LU Factorizations
We see that the G.E. algorithm can be used to actually calculate the LU factorization. We will look
at this in more detail in another example. We now examine how we can use the LU factorization
to solve the equation
Ax = b,
Since we have A = LU, we first solve
Iz = b,
then solve
\Jx = z.
Since L is lower triangular, we can solve for z with a forward substitution. Similarly, since U is
upper triangular, we can solve for x with a back substitution. We drag out the previous example
(which we never got around to solving):
/ 2
1
1
3
7 \
4
4
7
11
6
5
4
17
31
V 2
-1
7
15 J
We had found the LU factorization of A as
A =
So we solve
We get
Now we solve
We get the ugly solution
1
"
" 2
1
1
3
2
1
2
-2
1
3
1
1
3
7
1
-1
-1
1 _
_
12
" 1
7 "
2
1
11
3
1
1
z =
31
1
-1
-1
1
15 _
7
-3
2
1
1
3
2
-2
1
3
7
12
r 37 i
3.4. ITERATIVE SOLUTIONS
37
3.3.3 Some Theory
We aren't doing much proving here. The following theorem has an ugly proof in the Cheney &
Kincaid [7].
Theorem 3.3. If A is an n x n matrix, and naive Gaussian Elimination does not encounter a zero
pivot, then the algorithm generates a LU factorization of A, where L is the lower triangular part of
the output matrix, and U is the upper triangular part.
This theorem relies on us using the fancy version of G.E., which saves the multipliers in the
spots where there should be zeros. If correctly implemented, then, L is the lower triangular part
but with ones put on the diagonal.
This theorem is proved in Cheney & Kincaid [7]. This appears to me to be a case of something
which can be better illustrated with an example or two and some informal investigation. The proof
is an unillustrating index-chase-read it at your own risk.
3.3.4 Computing Inverses
We consider finding the inverse of A. Since
AA" 1 = I,
then the j th column of the inverse A -1 solves the equation
Ace — c j .
where ej is the column matrix of all zeros, but with a one in the j th position.
Thus we can find the inverse of A by running n linear solves. Obviously we are only going
to run G.E. once, to put the matrix in LU form, then run n solves using forward and backward
substitutions.
3.4 Iterative Solutions
Recall we are trying to solve
Ax = b.
We examine the computational cost of Gaussian Elimination to motivate the search for an alter-
native algorithm.
3.4.1 An Operation Count for Gaussian Elimination
We consider the number of floating point operations ("flops") required to solve the system Ate = b.
Gaussian Elimnation first uses row operations to transform the problem into an equivalent problem
of the form \Jx = b', where U is upper triangular. Then back substitution is used to solve for x.
First we look at how many floating point operations are required to reduce
/
an
ai2
ai 3 • •
O-ln
bi
\
a-21
a22
a23 • •
&2n
b 2
031
a33 • •
03n
bz
V
O-nl
a n 2
a«3 • •
b n
J
38
CHAPTER 3. SOLVING LINEAR SYSTEMS
to
Ol2
°22
a 32
ai3
a 23
a 33
a 2n
°3n
bi \
b's
\ a 'n2 «n3
i' 6' /
''nn u n /
First a multiplier is computed for each row. Then in each row the algorithm performs n
multiplies and n adds. This gives a total of (n— 1) + (n — l)n multiplies (counting in the computing
of the multiplier in each of the (n — 1) rows) and (n — l)n adds. In total this is 2n 2 — n — 1 floating
point operations to do a single pivot on the n by n system.
Then this has to be done recursively on the lower right subsystem, which is an (n— 1) by (n — 1)
system. This requires 2(n — 1) 2 — (n — 1) — 1 operations. Then this has to be done on the next
subsystem, requiring 2(n — 2) 2 — (n — 2) — 1 operations, and so on.
In total, then, we use /„ total floating point operations, with
n n n
j=l j=l j=l
Recalling that
n ^ n ^
^J 2 = ^W(n + l)(2n + l), and ^j = -(n)(n + l),
j'=i
We find that
1 2 o
/„ = -(4n - l)n(n + 1) - n w -n 6 .
b o
Now consider the costs of back substitution. To solve
/an • • • Ol,n-2 «l,n-l «ln \
V o
flra-2,n-2 O n _2,n-l On-2,7
a ra _i jn _i a n _i ; ,
a nn
bn-2
b n -i
bn j
for x n requires only a single division. Then to solve for x n -\ we compute
1
Xn—l
Q-n—l.n—l
[bn— 1 Q"n— l,n%n] >
and requires 3 flops. Similarly, solving for x n _2 requires 5 flops. Thus in total back substitution
requires B n total floating point operations with
B n = 2j — 1 = n(n — 1) — n = n(n — 2)mn z
3=1
3.4. ITERATIVE SOLUTIONS
39
3.4.2 Dividing by Multiplying
We saw that Gaussian Elimination requires around |ra 3 operations just to find the LU factorization,
then about n 2 operations to solve the system, when A is n x n. When n is large, this may take too
long to be practical. Additionally, if A is sparse (has few nonzero elements per row), we would like
the complexity of our computations to scale with the sparsity of A. Thus we look for an alternative
algorithm.
First we consider the simplest case, n = 1. Suppose we are to solve the equation
Ax = b.
for scalars A, b. We solve this by
\ 1 1
A ~ ujA ~ 1 - (1 - ojA)
uib =
1-r
uib,
where uo / is some real number chosen to "weight" the problem appropriately, and r = 1 — to A.
Now suppose that oj is chosen such that \r\ < 1. This can be done so long as A / 0, which would
have been a problem anyway. Now use the geometric expansion:
1 + r + r 2 + r 6 + . . .
Because of the assumption \r\ < 1, the terms r n converge to zero as n
approximate solution to our one dimensional problem as
oo. This gives the
1 + r + r 2 + r 3 + . . . + r k
Ldb +
Lob + r
r + r 2 +r 3 + ...+r k
ujb
ub
1 + r + r 2 +
+ r
k-l
cob
This suggests an iterative approach to solving Ax = b. First let x^ = cob, then let
x ( fc ) = ub + rx (fc_1) .
The iterates i^' will converge to the solution of Ax = b if \r\ < 1.
You should now convince yourself that because r n — > 0, that the choice of the initial iterate x^ ^
was immaterial, i.e., that under any choice of initial iterate convergence is guaranteed.
We now translate this scalar result into the vector case. The algorithm proceeds as follows: first
fix some initial estimate of the solution, x^°\ A good choice might be cob, but this is not necessary.
Then calculate successive approximations to the actual solution by updates of the form
x
W =U jb+ (I -ujk)x
(k-l)
It turns out that we can consider a slightly more general form of the algorithm, one in which
successive iterates are defined implicitly. That is we consider iterates of the form
Qa;( fc+1 ) = (Q - ljA) ajW + cob,
(3.3)
for some matrix Q, and some scaling factor to. Note that this update relies on vector additions
and possibly by premultiplication of a vector by A or Q. In the case where these two matrices are
sparse, such an update can be relatively cheap.
40
CHAPTER 3. SOLVING LINEAR SYSTEMS
Now suppose that as k — > oo, x^ converges to some vector x* , which is a fixed point of the
iteration. Then
Qx* = (Q -ujA)x* + ujb,
Qx* = Qx* - ujAx* + ujb,
u)Ax* = oob,
Ax* = b.
We have some freedom in choosing Q, but there are two considerations we should keep in mind:
1. Choice of Q affects convergence and speed of convergence of the method. In particular, we
want Q to be similar to A.
2. Choice of Q affects ease of computing the update. That is, given
z = (Q - A) + b,
we should pick Q such that the equation
Qa;(fc+i) = z
is easy to solve exactly.
These two goals conflict with each other. At one end of the spectrum is the so-called "impossible
iteration," at the other is the Richardsons.
3.4.3 Impossible Iteration
I made up the term "impossible iteration." But consider the method which takes Q to be A. This
seems to be the best choice for satisfying the first goal. Letting uj = 1, our method becomes
Aa ,(fc+i) = (A — A) £c( fe ) + b = b.
This method should clearly converge in one step. However, the second goal is totally ignored.
Indeed, we are considering iterative methods because we cannot easily solve this linear equation in
the first place.
3.4.4 Richardson Iteration
At the other end of the spectrum is the Richardson Iteration, which chooses Q to be the identity
matrix. Solving the system
Q^ +1 ) = z
is trivial: we just have x^ k+1 ^ = z.
Example Problem 3.4. Use Richardson Iteration with w = lon the system
" 6
1
1 "
" 12 "
A =
2
4
b =
1
2
6
6
Solution: We let
" l
"
" -5
-1
-1
Q =
1
, (Q-A) =
-2
-3
1
-1
-2
-5
3.4. ITERATIVE SOLUTIONS
41
We start with an arbitrary x^\ say = [2 2 2] T . We get x^ = [-2 - 10 - 10] T , and x^ =
[42 34 78] T .
Note the real solution is x = [2 — 1 1] T . The Richardson Iteration does not appear to converge
for this example, unfortunately. H
Example Problem 3.5. Apply Richardson Iteration with w = 1/6 on the previous system.
Solution: Our iteration becomes
-1/6
-1/6 "
' 2 "
x (k+l) =
-1/3
1/3
x^ +
_ -1/6
-1/3
1
We start with the same x^ as previously, x^ = [2 2 2] T . We get x^ = [4/3 0] T , x^ =
[2 - 4/9 7/9] T , and finally x^ = [2 - 0.99998 0.99998] T .
Thus, the choice of u has some affect on convergence. H
We can rethink the Richardson Iteration as
X ( k +V = (| _ W A) cc^ + cub = + w (b - AajW) .
Thus at each step we are adding some scaled version of the residual, defined as b — Ax , to the
iterate.
3.4.5 Jacobi Iteration
The Jacobi Iteration chooses Q to be the matrix consisting of the diagonal of A. This is more
similar to A than the identity matrix, but nearly as simple to invert.
Example Problem 3.6. Use Jacobi Iteration, with u> = 1, to solve the system
" 6
1
1 "
" 12 "
A =
2
4
1
b
1
2
6
6
We let
" 6
"
-1
-1 "
- 1
6
Q =
4
, (Q-
-A)
2
Q- l =
1
4
6
1
-2
We start with an arbitrary say x^ = [2 2 2]
[t ~ I T ] T • Continuing, we find that x^ w [1.987 - 1.019 0.981]
Note the real solution is x = [2 — 1 1] T .
We get x^ = [§ - 1 0] T . Then x^ =
T
There is an alternative way to describe the Jacobi Iteration for uj = 1. By considering the update
elementwise, we see that the operation can be described by
x
(fe+i) _ 1
Thus an update takes less than 2n 2 operations. In fact, if A is sparse, with less than k nonzero
entries per row, the update should take less than 2nk operations.
42
CHAPTER 3. SOLVING LINEAR SYSTEMS
3.4.6 Gauss Seidel Iteration
The Gauss Seidel Iteration chooses Q to be lower triangular part of A, including the diagonal. In
this case solving the system
Qx^ = z
is performed by forward substitution. Here the Q is more like A than for Jacobi Iteration, but
involves more work for inverting.
Example Problem 3.7. Use Gauss Seidel Iteration to again solve for
" 6
1
1 "
" 12 "
A =
2
4
b =
1
2
6
6
Solution: We let
Q
6
2 4
1 2 6
(Q-A)
-1
-1
We start with an arbitrary x(°), say x<W = [2 2 2]. We get xW = [§ - § l] . Then x
[35 _ 35 il T
L 18 36 1 1 ■
Already this is fairly close to the actual solution x = [2 — 1 1] T .
,(2)
Just as with Jacobi Iteration, there is an easier way to describe the Gauss Seidel Iteration. In
this case we will keep a single vector x and overwrite it, element by element. Thus for j = 1, 2, . . . , n,
we set
Xj *~ I °j 2s a i* x * I •
" \ i=l&j j
This looks exactly like the Jacobi update. However, in the sum on the right there are some "old"
values of Xj and some "new" values; the new values are those X{ for which i < j.
Again this takes less than 2n 2 operations. Or less than Ink if A is sufficiently sparse.
An alteration of the Gauss Seidel Iteration is to make successive "sweeps" of this redefinition,
one for j = 1, 2, . . . , n, the next for j = n, n — 1, . . . , 2, 1. This amounts to running Gauss Seidel
once with Q the lower triangular part of A, then running it with Q the upper triangular part. This
iterative method is known as "red-black Gauss Seidel."
3.4.7 Error Analysis
Suppose that x is the solution to equation 3.4. Define the error vector:
3.4. ITERATIVE SOLUTIONS
43
Now notice that
x
x
X
(fc+1) =
(fc+1) =
(fc+1) _
Q _1 (Q - ujA) £c (fe) + Q _1 u;b,
Q-^ajW - wQ _1 AajW + wQ^Aa;,
a?
(fc+i)
a; = a;
(fc)
IE
(fc)
a;
- x - wQ^A (x^ - x^j ,
(fc)
= (l-wQ^A) cW.
e (fe+i)
Reusing this relation we find that
e« = (I - wQ^A) e^" 1 ),
= (l-cQ- 1 A) fc e(°).
We want to ensure that e^ k+1 ^ is "smaller" than e^ k \ To do this we recall matrix and vector norms
from Subsection 1.4.1.
,(*)
(l-wQ- 1 A) fc c(°) < j|l -t^Q^A
,(o)
(See Example Problem 1.29.)
Thus our iteration converges (e^ goes to the zero vector, i.e., x^ — ► x) if
|l -wQ _1 A|| 2 < 1.
This gives the theorem:
Theorem 3.8. An iterative solution scheme converges for any starting x^ if and only if all
eigenvalues of I — wQ _1 A are less than 1 in absolute value, i.e., if and only if
||l -<jQ _1 A|| 2 < 1
Another way of saying this is "the spectral radius of I — wQ _1 A is less than 1."
In fact, the speed of convergence is decided by the spectral radius of the matrix-convergence
is faster for smaller values. Recall our introduction to iterative methods in the scalar case, where
the result relied on uj being chosen such that |1 — ujA\ < 1. You should now think about how
eigenvalues generalize the absolute value of a scalar, and how this relates to the norm of matrices.
Let y be an eigenvector for Q _1 A, with corresponding eigenvalue A. Then
(I - wQ _1 A) y = y- ujQ- l Ay = y-uXy = (l- uX) y.
This relation may allow us to pick the optimal uj for given A,Q. It can also show us that
sometimes no choice of uj will give convergence of the method. There are a number of different
related results that show when various methods will work for certain choices of uj. We leave these
to the exercises.
44
CHAPTER 3. SOLVING LINEAR SYSTEMS
Example Problem 3.9. Find conditions on uj which guarantee convergence of Richardson's Iter-
ation for finding approximate iterative solutions to the system Ax = b, where
" 6
1
1 "
" 12 "
A =
2
4
b =
1
2
6
6
Solution: By Theorem 3.8, with Q the identity matrix, we have convergence if and only if
||l - wA|| 2 < 1
We now use the fact that "eigenvalues commute with polynomials;" that is if f(x) is a polynomial
and A is an eigenvalue of a matrix A, then /(A) is an eigenvalue of the matrix /(A). In this case the
polynomial we consider is f(x) = x°—lux 1 . Using octave or Matlab you will find that the eigenvalues
of A are approximately 7.7321,4.2679, and 4. Thus the eigenvalues of I — uiA are approximately
1 - 7.7321w, 1 - 4.2679cj, 1 - 4u.
With some work it can be shown that all three of these values will be less than one in absolute
value if and only if
< uj < — ^ — « 0.388
7.7321
(See also Exercise (3.10).)
Compare this to the results of Example Problem 3.4, where for this system, ui = 1 apparently did
not lead to convergence, while for Example Problem 3.5, with uj = 1/6, convergence was observed.
H
3.4.8 A Free Lunch?
The analysis leading to Theorem 3.8 leads to an interesting possible variant of the iterative scheme.
For simplicity we will only consider an alteration of Richardson's Iteration. In the altered algorithm
we presuppose the existence, via some oracle, of a sequence of weightings, u>i, which we use in each
iterative update. Thus our algorithm becomes:
1. Select some initial iterate x^°\
2. Given iterate £c( fc_1 ), define
X W = (\ - uj k A) x^ + u; k b.
Following the analysis for Theorem 3.8, it can be shown that
e ( fc ) = (I - u k A) e^ -1 )
where, again, = x^ — x, with x the actual solution to the linear system. Expanding e^ -1 )
similarly gives
e W = ( |_^ A )(l^_ 1 A)e^ 2 )
3.4. ITERATIVE SOLUTIONS
45
Remember that we want to be small in magnitude, or, better yet, to be zero. One way to
guarantee that is zero, regardless of the choice of e(°) it to somehow ensure that the matrix
k
B = H I - u>iA
1=1
has all zero eigenvalues.
We again make use of the fact that eigenvalues "commute" with polynomials to claim that if
Xj is an eigenvalue of A, then
k
n 1 - ^
i=i
is an eigenvalue of B. This eigenvalue is zero if one of the uj{ for 1 < % < k is 1/Aj. This suggests
how we are to pick the weightings: let them be the inverses of the eigenvalues of A.
In fact, if A has a small number of distinct eigenvalues, say m eigenvalues, then convergence to
the exact solution could be guaranteed after only m iterations, regardless of the size of the matrix.
As you may have guessed from the title of this subsection, this is not exactly a practical
algorithm. The problem is that it is not simple to find, for a given arbitrary matrix A, one, some,
or all its eigenvalues. This problem is of sufficient complexity to outweigh any savings to be gotten
from our "free lunch" algorithm.
However, in some limited situations this algorithm might be practical if the eigenvalues of A
are known a priori.
46
CHAPTER 3. SOLVING LINEAR SYSTEMS
Exercises
(3.1) Find the LU decomposition of the following matrices, using naive Gaussian Elimination
(a)
3 -1 -2
9 -1 -4
-6 10 13
(b)
24
12
13
16
11
19
-3
1
-4
6
-2
5 -8
(3.2) Perform back substitution to solve the equation
1
3
5
3 "
' -1 "
2
4
3
-1
x =
2
1
1
2
2
(3.3) Perform Naive Gaussian Elimination to prove Cramer's rule for the 2D case. That is, prove
that the solution to
a b
X
" / "
c d
. v _
. a _
is given by
det
det
a f
c g
a b
c d
and x
det
det
/ b
9 d
a b
c d
(3.4)
(3.5)
(3.6)
Implement Cramer's rule to solve a pair of linear equations in 2 variables. Your m-file should
have header line like:
function x = cramer2(A,b)
where A is a 2 x 2 matrix, and b and x are 2x1 vectors. Your code should find the x such
that Ax = b. (See Exercise (3.3))
Test your code on the following (augmented) systems:
(a)
(d)
3 -2
4 -3
-0.0590
0.1080
(b)
0.2372
-0.4348
1.24
-0.744
-0.3528
0.6452
-3.48
2.088
(c)
1.24 -3.48
-0.744 2.088
1
-0.6
Given two lines parametrized by f(t) = at + b, and g(s) = cs + d, set up a linear 2x2 system
of equations to find the t, s at the point of intersection of the two lines. If you were going
to write a program to detect the intersection of two lines, how would you detect whether
they are parallel? How is this related to the form of the solution to a system of two linear
equations? (See Exercise (3.3))
Prove that the inverse of the matrix
1 •••
02 1 ■■■
3 1 •■■
a n ■■■ 1
3.4. ITERATIVE SOLUTIONS
47
is the matrix
1
-a 2 1
-a 3 1
-Or,
1
(Hint: Multiply them together.)
(3.7) Under the strategy of scaled partial pivoting, which row of the following matrix will be the
first pivot row?
10
17
-10
0.1
0.9
-3
3
-3
0.3
-4
0.3
0.1
0.01
-1
0.5
2
3
4
-3
5
10
100
1
0.1
(3.8) Let A be a symmetric positive definite n x n matrix with n distinct eigenvalues. Letting
y(°) = 5/ ||6|| 2 5 consider the iteration
y
(fc+i) _
WW
(a) What is ||y (fc) || 2 ?
(b) Show that y( fc )
(c) Show that as k
A k b/\\A k b\\ 2 .
oo, y^> converges to the (normalized) eigenvector associated with the
largest eigenvalue of A.
(3.9) Consider the equation
Letting x^ = [11 0] T , find the iterate a^ 1 ) by one step of Richardson's Method. And by
one step of Jacobi Iteration. And by Gauss Seidel.
(3.10) Let A be a symmetric n x n matrix with eigenvalues in the interval [a, with < a < (3,
and a + (3 ^ 0. Consider Richardson's Iteration
1
3
5
' -5 "
-2
2
4
x =
-6
4
-3
-4
10
X
(fc+1)
(I — uA) x^> + ub.
Recall that e( fc+1 ) = (I - ujA) .
(a) Show that the eigenvalues of I — ujA are in the interval [1 — u>(3, 1 — wa].
(b) Prove that
max {|A| : 1 — tu(3 < A < 1 — iva}
ujI3
is minimized when we choose ui such that 1
look at the graph of something versus to.)
(c) Show that this relationship is satisfied by u> = 2/ (a + 0).
(d) For this choice of u show that the spectral radius of I — ooA is
\a + [3\'
(1 — ua) . (Hint: It may help to
48
CHAPTER 3. SOLVING LINEAR SYSTEMS
(3.11
(e) Show that when < a, this quantity is always smaller than 1.
(f) Prove that if A is positive definite, then there is an uj such that Richardson's Iteration
with this lo will converge for any choice of x (°) .
(g) For which matrix do you expect faster convergence of Richardson's Iteration: Ai with
eigenvalues in [10,20] or A2 with eigenvalues in [1010, 1020]? Why?
I Implement Richardson's Iteration to solve the system Ax = b. Your m-file should have
header line like:
function xk = richardsons(A,b,xO,w,k)
Your code should return x^ based on the iteration
x
O'+i)
X
(i)
UJ
(Ax®
Let w take the place of uj, and let xO be the initial iterate x^°\ Test your code for A, b for
which you know the actual solution to the problem. (Hint: Start with A and the solution x
and generate b.) Test your code on the following matrices:
• Let A be the Hilbert Matrix. This is generated by the octave command A = hilb(10),
or whatever (smallish) integer. Try different values of uj, including uj = 1.
• Let A be a Toeplitz matrix of the form:
A
These can be generated by the octave command A = toeplitz ([-2 100000 0]).
Try different values of uj, including uj = —1/2.
(3.12) Let A be a nonsingular n x n matrix. We wish to solve Ax = b. Let x^ be some starting
vector, let Vk be span{r^,Ar^, . . . , A fc r^} , and let Vk be the set of polynomials, p(x) of
degree k with p(0) = 1.
Consider the following iterative method: Let cc( fe+1 ) be the x that solves
-2
1
••
•
1
-2
1 ••
•
1
-2 ••
•
••
• -2
mm
\b - Ax\
2 •
Let rW = b- Ax^ k \
(a) Show that if x € a;( ) + V k , then b - Ax = p(A)r^ for some p G V k .
(b) Prove that, conversely, for any p G Vk there is some x G x(°) + Vk, such that b
p(A)r(°).
(c) Argue that
Ax
,(fc+i)
mm
pev k
p(A)r
(0)
(d) Prove that this iteration converges in at most n steps. (Hint: Argue for the existence
of a polynomial in V n that vanishes at all the eigenvalues of A. Use this polynomial to
show that ||r( n )|L < 0.)
Chapter 4
Finding Roots
4.1 Bisection
We are now concerned with the problem of finding a zero for the function f(x), i.e., some c such
that /(c) = 0.
The simplest method is that of bisection. The following theorem, from calculus class, insures
the success of the method
Theorem 4.1 (Intermediate Value Theorem). If f(x) is continuous on [a,b] then for any y
such that y is between f(a) and f(b) there is a c G [a, b] such that f(c) = y.
The IVT is best illustrated graphically. Note that continuity is really a requirement here-a
single point of discontinuity could ruin your whole day, as the following example illustrates.
Example 4.2. The function f(x) = \ is not continuous at 0. Thus if G [a, b] , we cannot apply
the IVT. In particular, if G [a,b] it happens to be the case that for every y between f(a),f(b)
there is no c G [a, b] such that /(c) = y.
In particular, the IVT tells us that if f(x) is continuous and we know a, b such that /(a), f(b)
have different sign, then there is some root in [a,b]. A decent estimate of the root is c =
We can check whether /(c) = 0. If this does not hold then one and only one of the two following
options holds:
1. /(a), /(c) have different signs.
2. /(c), f(b) have different signs.
We now choose to recursively apply bisection to either [a,c] or [c, b], respectively, depending on
which of these two options hold.
4.1.1 Modifications
Unfortunately, it is impossible for a computer to test whether a given black box function is contin-
uous. Thus malicious or incompetent users could cause a naively implemented bisection algorithm
to fail. There are a number of easily conceivable problems:
1. The user might give f,a,b such that f(a),f(b) have the same sign. In this case the function
/ might be legitimately continuous, and might have a root in the interval [a, b] . If, taking
c = /(a), f(b), /(c) all have the same sign, the algorithm would be at an impasse. We
should perform a "sanity check" on the input to make sure f(a),f(b) have different signs.
49
50
CHAPTER 4. FINDING ROOTS
2. The user might give f,a,b such that / is not continuous on [a, b], moreover has no root in
the interval [a, b] . For a poorly implemented algorithm, this might lead to an infinite search
on smaller and smaller intervals about some discontinuity of /. In fact, the algorithm might
descend to intervals as small as machine precision, in which case the midpoint of the interval
will, due to rounding, be the same as one of the endpoints, resulting in an infinite recursion.
3. The user might give / such that / has no root c that is representable in the computer's
memory. Recall that we think of computers as storing numbers in the form ±r x 10 fc ; given
a finite number of bits to represent a number, only a finite number of such numbers can be
represented. It may legitimately be the case that none of them is a root to /. In this case,
the behaviour of the algorithm may be like that of the previous case. A well implemented
version of bisection should check the length of its input interval, and give up if the length is
too small, compared to machine precision.
Another common error occurs in the testing of the signs of /(a), f(b). A slick programmer might
try to implement this test in the pseudocode:
if (f (a)f (b) > 0) then . . .
Note however, that |/(a)| , |/(&)| might be very small, and that f(a)f(b) might be too small to
be representable in the computer; this calculation would be rounded to zero, and unpredictable
behaviour would ensue. A wiser choice is
if (sign(f(a)) * sign(f(b)) > 0) then ...
where the function sign (x) returns —1,0, 1 depending on whether x is negative, zero, or positive,
respectively.
The pseudocode bisection algorithm is given as Algorithm 1.
4.1.2 Convergence
We can see that each time recursive_bisection (/, a,b, . . .) is called that \b — a\ is half the length
of the interval in the previous call. Formally call the first interval [ao,&oL an d the first midpoint
Co- Let the second interval be [ai,&i], etc. Note that one of a±,bi will be cp, and the other will be
either ao or bo. We are claiming that
, _ bn-i — a n -i
bp - ap
~ 2 n
Theorem 4.3 (Bisection Method Theorem). If f(x) is a continuous function on [a, b] such that
f(a)f(b) < 0, then after n steps, the algorithm run_bisection will return c such that \c — c'\ <
^ft^, where d is some root of /.
4.2 Newton's Method
Newton's method is an iterative method for root finding. That is, starting from some guess at the
root, xp, one iteration of the algorithm produces a number x\, which is supposed to be closer to a
root; guesses X2,x%, . . . ,x n follow identically.
Newton's method uses "linearization" to find an approximate root. Recalling Taylor's Theorem,
we know that
f(x + h)^f(x) + f'(x)h.
4.2. NEWTON'S METHOD
51
Algorithm 1: Algorithm for finding root by bisection.
Input: a function, two endpoints, a x-tolerance, and a y-tolerance
Output: a c such that \f(c)\ is smaller than the y-tolerance.
run_bisection(/, a, b, 5, e)
(1) Let /a <- sign (/(a)).
(2) Let fb<- sign (/(&)).
(3) if /a/6 >
(4) throw an error.
(5) else if fafb =
(6) if fa =
(7) return a
(8) else
(9) return b
(10) Let c ^ ^
(11) while 6-a>
(12) Let /c <- /(c).
(13) if |/c| < e
(14) return c
(15) if sign (fa) sign (/c) <
(16) Let b^c, /&<-/c, c ^- s±£.
(17) else
(18) Let a <- c, fa <- /c, c «-
(19) return c
This approximation is better when /"(•) is "well-behaved" between x and x + h. Newton's method
attempts to find some h such that
This is easily solved as
= f(x + h) = f(x) + f'(x)h.
-fix)
h
fix)
An iteration of Newton's method, then, takes some guess Xk and returns Xk+i defined by
Xk+l = x k
f(Xk)
f'(x k y
(4.1)
An iteration of Newton's method is shown in Figure 4.1, along with the linearization of f(x) at
Xk-
4.2.1 Implementation
Use of Newton's method requires that the function fix) be differentiable. Moreover, the derivative
of the function must be known. This may preclude Newton's method from being used when f(x)
is a black box. As is the case for the bisection method, our algorithm cannot explicitly check for
continuity of f(x). Moreover, the success of Newton's method is dependent on the initial guess
xq. This was also the case with bisection, but for bisection there was an easy test of the initial
interval-i.e., test if /(ao)/(&o) < 0-
52
CHAPTER 4. FINDING ROOTS
Figure 4.1: One iteration of Newton's method is shown for a quadratic function f(x). The lin-
earization of f{x) at Xk is shown. It is clear that xt+i is a root of the linearization. It happens to
be the case that \f(xk+i)\ is smaller than , i.e., x^+i is a better guess than xp..
Our algorithm will test for goodness of the estimate by looking at • The algorithm will
also test for near-zero derivative. Note that if it were the case that f'( x k) = then h would be ill
defined.
Algorithm 2: Algorithm for finding root by Newton's Method.
Input: a function, its derivative, an initial guess, an iteration limit, and a tolerance
Output: a point for which the function has small value.
run_newton(/, /', xO, N, tol)
(1)
Let x <— xO, n <— 0.
(2)
while n < N
(3)
Let fx <- f{x).
(4)
if \fx\ < tol
(5)
return x.
(6)
Let fpx <— f'(x).
(7)
if \fpx\ < tol
(8)
Warn u f'(x) is small
(9)
return x.
(10)
Let x <— x — fx/ fpx.
(11)
Let n <— n + 1.
4.2. NEWTON'S METHOD
53
4.2.2 Problems
As mentioned above, convergence is dependent on /(x), and the initial estimate xq. A number of
conceivable problems might come up. We illustrate them here.
Example 4.4. Consider Newton's method applied to the function /(x) = x-? with j > 1, and with
initial estimate Xq ^ 0.
Note that f(x) = has the single root x = 0. Now note that
4 ^ 1
Xk+l =X k 73T = 1 - - )x k .
jx J k V 3 J
Since the equation has the single root x = 0, we find that x& is converging to the root. However,
it is converging at a rate slower than we expect from Newton's method: at each step we have a
constant decrease of 1 — (l/j) , which is a larger number (and thus worse decrease) when j is larger.
Example 4.5. Consider Newton's method applied to the function f(x) = with initial estimate
x = 3.
Note that /(x) is continuous on R + . It has a single root at x = 1. Our initial guess is not too far
from this root. However, consider the derivative:
,,, . x\ - lnx 1 - lnx
fix) = ~^— 2 = —
x z x z
If x > e 1 , then 1 — lnx < 0, and so /'(x) < 0. However, for x > 1, we know /(x) > 0. Thus taking
Xk+l =X k - -J7-, r > X k .
f'(Xk)
The estimates will "run away" from the root x = 1.
Example 4.6. Consider Newton's method applied to the function f(x) = sin (x) for the initial
estimate xo 7^ 0, where xq has the odious property 2xo = tanxo-
You should verify that there are an infinite number of such xq. Consider the identity of X\ :
f(x ) sin(x )
Xl = X - — = Xq = Xq- tanx = X - 2x = -X .
f'(x ) cos(xo)
Now consider X2 :
/(xi) sin(— xo) sin(xo)
x 2 = xi - — — - = -x 7 r - x H 7 — r = -x + tan x = -x + 2x = £o-
/'(xi) cos(-xo) cos(x )
Thus Newton's method "cycles" between the two values xo, — xo-
Of course, Newton's method may find some iterate Xk for which f'(x k ) = 0, in which case, there
is no well-defined x^+i.
4.2.3 Convergence
When Newton's Method converges, it actually displays quadratic convergence. That is, if e k = x^— r,
1 1 2
where r is the root that the Xk are converging to, that |efc+i| < C | | . If, for example, it were
the case that C = 1, then we would double the accuracy of our root estimate with each iterate.
That is, if eo were 0.001, we would expect e\ to be on the order of 0.000001. The following theorem
formalizes our claim:
54
CHAPTER 4. FINDING ROOTS
Theorem 4.7 (Newton's Method Convergence). If f(x) has two continuous derivatives, and
r is a simple root of f(x), then there is some D such that if \xq — r\ < D, Newton's method will
converge quadratically to r.
The proof is based on arguments from real analysis, and is omitted; see Cheney &; Kincaid for
the proof [7]. Take note of the following, though:
1. The proof requires that r be a simple root, that is that f'(r) ^ 0. When this does not hold
we may get only linear convergence, as in Example 4.4.
2. The key to the proof is using Taylor's theorem, and the definition of Newton's method to
show that
k+1 2f'(x k ) '
where is between Xk and r = Xk + e&.
The proof then proceeds by showing that there is some region about r such that
(a) in this region \ f"{x)\ is not too large and is not too small, and
(b) if Xk is in the region, then Xk+i is in the region.
In particular the region is chosen to be so small that if Xk is in the region, then the factor
e\ will outweigh the factor /^\f'( x k)\ ■ You can roughly think of this region as an
"attractor" region for the root.
3. The theorem never guarantees that some Xk will fall into the "attractor" region of a root r of
the function, as in Example 4.5 and Example 4.6. The theorem that follows gives sufficient
conditions for convergence to a particular root.
Theorem 4.8 (Newton's Method Convergence II [2]). If f(x) has two continuous derivatives
on [a,b], and
1. f(a)f(b) < 0,
2. /'(x)^0on [a, b],
3. fix) does not change sign on [a, b],
4. Both |/(a)| < (6 - a) |/'(a)| and |/(6)| < (b - a) \f'(b)\ hold,
then Newton's method converges to the unique root of f(x) = for any choice of xq £ [a, b] .
4.2.4 Using Newton's Method
Newton's Method can be used to program more complex functions using only simple functions.
Suppose you had a computer which could perform addition, subtraction, multiplication, division,
and storing and retrieving numbers, and it was your task to write a subroutine to compute some
complex function g(-). One way of solving this problem is to have your subroutine use Newton's
Method to solve some equation equivalent to g(z) — x = 0, where z is the input to the subroutine.
Note that it is assumed the subroutine cannot evaluate g(z) directly, so this equation needs to be
modified to satisfy the computer's constraints.
Quite often when dealing with problems of this type, students make the mistake of using New-
ton's Method to try to solve a linear equation. This should be an indication that a mistake was
made, since Newton's Method can solve a linear equation in a single step:
Example 4.9. Consider Newton's method applied to the function f(x) = ax + b. The iteration is
given as
axk + b
x k+ i <- x k •
a
This can be rewritten simply as Xk+i < b/a.
4.3. SECANT METHOD
55
The following example problems should illustrate this process of "bootstrapping" via Newton's
Method.
Example Problem 4.10. Devise a subroutine using only subtraction and multiplication that can
find the multiplicative inverse of an input number z, i.e., can find (1/z) .
Solution: We are tempted to use the linear function f(x) = (1/z) — x. But this is a linear equation
for which Newton's Method would reduce to Xk+i (1/z) • Since the subroutine can only use
subtraction and multiplication, this will not work.
Instead apply Newton's Method to f(x) = z — (1/x) . The Newton step is
z — (1/xk) 2
x k+l ^~ x k 7Z i 2\ = x k ~ ZX k + Xk = x k (2 — ZXk) ■
Note this step uses only multiplication and subtraction. The subroutine is given in Algorithm 3.
H
Algorithm 3: Algorithm for finding a multiplicative inverse using simple operations.
Input: a number
Output: its multiplicative inverse
INVS(z)
(1) if z = throw an error.
(2) Let x <— sign (z) , n <— 0.
(3) while n < 50
(4) Let x <- x (2 - zx) .
(5) Let n <— n + 1.
Example Problem 4.11. Devise a subroutine using only simple operations which computes the
square root of an input number z.
Solution: The temptation is to find a zero for f(x) = \fz — x. However, this equation is linear
in x. Instead let f(x) = z — x 2 . You can easily see that if x is a positive root of f(x) = 0, then
x = \fz. The Newton step becomes
x k+ i <- x k
4
-2x k
after some simplification this becomes
x k , z
Xk+1 ^ T 2x~ k
Note this relies only on addition, multiplication and division.
The final subroutine is given in Algorithm 4.
H
4.3 Secant Method
The secant method for root finding is roughly based on Newton's method; however, it is not
assumed that the derivative of the function is known, rather the derivative is approximated by
the value of the function at some of the iterates, x k - More specifically, the slope of the tangent
56
CHAPTER 4. FINDING ROOTS
Algorithm 4: Algorithm for finding a square root using simple operations.
Input: a number
Output: its square root
SQRt(z)
(1) if z < throw an error.
(2) Let x <- 1, n <- 0.
(3) while n < 50
(4) Let x <- (a; + z/x) /2.
(5) Let n n + 1.
Figure 4.2: One iteration of the Secant method is shown for some quadratic function f(x). The
secant line through (x k -i, f(xk~i)) and (xk,f(xk)) is shown. It happens to be the case that
is smaller than \f(xk)\ , i.e., is a better guess than x^.
line at (x^, f(xk)) , which is f'(xk) is approximated by the slope of the secant line passing through
(x k -i,f(xk-i)) and (x k ,f(x k )) , which is
f(Xk) ~ f(Xk-l)
%k %k— 1
Thus the iterate x/c+i is the root of this secant line. That is, it is a root to the equation
— {x-x k ) = y- f{x k ).
4.3. SECANT METHOD
Since the root has a y value of 0, we have
57
fjxk) ~ fi x k-
Xk Xk—l
(x k+1 - x k ) = -f(x k ),
_ _ ( x k - Xk-1
\f(xk) - f(x k -
1,
x k+ i = x k -
Xk X k —i
f(xk) - f(Xk-l)
f(x k ),
f{x k ).
(4.2)
You will note this is the recurrence of Newton's method, but with the slope f'(x k ) substituted
by the slope of the secant line. Note also that the secant method requires two initial guesses, xo, x±,
but can be used on a black box function. The secant method can suffer from some of the same
problems that Newton's method does, as we will see.
An iteration of the secant method is shown in Figure 4.2, along with the secant line.
Example 4.12. Consider the secant method used on x 3 + x 2 — x — 1, with xo = 2, x\
Note that this function is continuous and has roots ±1. We give the iterates here:
A:
Xk
f(xk)
2
9
1
0.5
-1.125
2
0.666666666666667
-0.925925925925926
3
1.44186046511628
2.63467367653163
4
0.868254072087394
-0.459842466254495
5
0.953491494113659
-0.177482458876898
6
1.00706900811804
0.0284762692197613
7
0.999661272951803
-0.0013544492875992
8
0.999997617569723
-9.52969840528617 x 10 - '
9
1.0000000008072
3.22880033820638 x 10~ 09
10
0.999999999999998
-7.43849426498855 x 10 _
4.3.1 Problems
As with Newton's method, convergence is dependent on f(x), and the initial estimates xq,X\. We
illustrate a few possible problems:
Example 4.13. Consider the secant method for the function f(x) = with x$ = 3, X\ = 4.
As with Newton's method, the iterates diverge towards infinity, looking for a nonexistant root. We
58
CHAPTER 4. FINDING ROOTS
rive some iterates here:
k Xfc
3
1 4
2 21.6548475770851
3 33.9111765137635
4 67.3380435135758
5 117.820919458675
6 210.543986613847
7 366.889164762149
8 637.060241341843
9 1096.54125113444
10 1878.34688714646
11 3201.94672271613
12 5437.69020766155
13 9203.60222260594
14 15533.1606791089
15 26149.7196085218
16 43924.8466075548
17 73636.673898472
f(xk)
0.366204096222703
0.346573590279973
0.142011128224341
0.103911011441661
0.0625163004418104
0.0404780904944712
0.0254089165873003
0.0160949419231219
0.010135406045582
0.00638363233543847
0.00401318169994875
0.0025208146648422
0.00158175793894727
0.000991714984152597
0.000621298692241343
0.000388975250950428
0.000243375589137882
0.000152191807070607
Example 4.14. Consider Newton's method applied to the function f(x) = — — jj, with initial
estimates xq = — l,xi = 1.
You can easily verify that /(#o) = f(xi), and thus the secant line is horizontal. And thus X2 is not
defined.
Algorithm 5: Algorithm for finding root by secant method.
Input: a function, initial guesses, an iteration limit, and a tolerance
Output: a point for which the function has small value.
run_secant(/, xO, xl, N, tol)
(1) Let x <— xl, xp <— xO, fh <— /(xl), fp <— f(x0),n <— 0.
(2) while n < N
(3) if \fh\ < tol
(4) return x.
(5) Let fpx <- (fh - fp)/ (x - xp).
(6) if \fpx\< tol
(7) Warn "secant slope is too small; giving up."
(8) return x.
(9) Let xp <— x, fp <— //t.
(10) Let x <— x — fh/ fpx.
(11) Let fh<-f(x). '
(12) Letnf-n + 1.
4.3.2 Convergence
We consider convergence analysis as in Newton's method. We assume that r is a root of f(x), and
let e ra = r — x n . Because the secant method involves two iterates, we assume that we will find some
4.3. SECANT METHOD
59
relation between e k+ \ and the previous two errors, e k ,e k -i.
Indeed, this is the case. Omitting all the nasty details (see Cheney & Kincaid [7]), we arrive at
the imprecise equation:
2/'(r)
Again, the proof relies on finding some "attractor" region and using continuity.
We now postulate that the error terms for the secant method follow some power law of the
following type:
|e fc+ i| ~ A\e k \ a .
Recall that this held true for Newton's method, with a = 2. We try to find the a for the secant
method. Note that
II A I \C
\ek\ = A\e k -i\ ,
so
I I A-± I |i
|efe_i| = A <* \e k \ a ■
Then we have
A \e k \ a = |e fc+ i| = C \e k \ |e fc _i| = CA~* |e fe | 1+ ~ ,
Since this equation is to hold for all e k , we must have
1
a = 1 + -.
a
This is solved by a = \ (l + V5) w 1.62. Thus we say that the secant method enjoys superlinear
convergence; This is somewhere between the convergence rates for bisection and Newton's method.
60
CHAPTER 4. FINDING ROOTS
Exercises
(4.1) Consider the bisection method applied to find the zero of the function f(x) = x 2 — 5x + 3,
with ao = 0, &o = 1- What are a\, b\l What are 02, 62 ?
(4.2) Approximate \/l0 by using two steps of Newton's method, with an initial estimate of xq = 3.
(c/. Example Problem 4.11) Your answer should be correct to 5 decimal places.
(4.3) Consider bisection for finding the root to cosx = 0. Let the initial interval Iq be [0, 2]. What
is the next interval considered, call it I\! What is I-p. I%1
(4.4) What does the sequence defined by
1 1
x = 1, Xfc+i = -x k H
2 x k
converge to?
(4.5) Devise a subroutine using only simple operations that finds, via Newton's Method, the cubed
root of some input number z.
(4.6) Use Newton's Method to approximate \^9. Start with xq = 2. Find X2.
(4.7) Use Newton's Method to devise a sequence xq,x\, . . . such that Xk — > In 10. Is this a reasonable
way to write a subroutine that, given z, computes In z? (Hint: such a subroutine would require
computation of e Xk . Is this possible for rational Xk without using a logarithm? Is it practical?)
(4.8) Give an example (graphical or analytic) of a function, f(x) for which Newton's Method:
(a) Does not find a root for some choices of xq.
(b) Finds a root for every choice of xq.
(c) Falls into a cycle for some choice of xq.
(d) Converges slowly to the zero of f(x).
(4.9) How will Newton's Method perform for finding the root to f(x) = y/\x\ = 0?
(4.10) Implement the inverse finder described in Example Problem 4.10. Your m-file should have
header line like:
function c = invs(z)
where z is the number to be inverted. You may wish to use the builtin function sign. As
an extra termination condition, you should have the subroutine return the current iterate if
zx is sufficiently close to 1, say within the interval (0.9999,0.0001). Can you find some z for
which the algorithm performs poorly?
(4.11) Implement the bisection method in octave/Matlab. Your m-file should have header line like:
function c = run_bisection(f , a, b, tol)
where f is the name of the function. Recall that feval(f ,x) when f is a string with the
name of some function evaluates that function at x. This works for builtin functions and
m-files.
(a) Run your code on the function f(x) = cosx, with ao = 0, 60 = 2. In this case you can
set f = "cos".
(b) Run your code on the function f(x) = x 2 — 5x + 3, with ao = 0, 60 = 1- In this case you
will have to set f to the name of an m-file (without the ".m") which will evaluate the
given f(x).
(c) Run your code on the function f(x) = x — cosx, with ao = 0, 60 = 1-
(4.12) Implement Newton's Method. Your m-file should have header line like:
function x = run_newton(f , fp, xO, N, tol)
where f is the name of the function, and f p is its derivative. Run your code to find zeroes of
the following functions:
(a) f(x) = tanx — x.
4.3. SECANT METHOD
61
(b) f(x) = x 2 — (2 + e)x + 1 + e, for e small.
(c) f(x) = x 7 - 7x 6 + 21a; 5 - 35x 4 + 35x 3 - 21x 2 + 7x-l.
(d) f(x) = (x-l) 7 .
(4.13) Implement the secant method Your m-file should have header line like:
function x = run_secant (f , xO, xl, N, tol)
where f is the name of the function. Run your code to find zeroes of the following functions:
(a) f(x)=t&nx — x.
(b) f{x) = x 2 + 1. (Note: This function has no zeroes.)
(c) f(x) = 2 + sinx. (Note: This function also has no zeroes.)
CHAPTER 4. FINDING ROOTS
Chapter 5
Interpolation
5.1 Polynomial Interpolation
We consider the problem of finding a polynomial that interpolates a given set of values:
X
x
Xn
y
yo
yi
Vn
where the Xi are all distinct. A polynomial p(x) is said to interpolate these data if p(xi) = yi for
i = 0, 1, . . . , n. The x« values are called "nodes."
Sometimes, we will consider a variant of this problem: we have some black box function, f(x),
which we want to approximate with a polynomial p{x). We do this by finding the polynomial
interpolant to the data
X
x
Xl
Xn
m
f(xo)
f(Xn)
for some choice of distinct nodes Xi.
5.1.1 Lagranges Method
As you might have guessed, for any such set of data, there is an n-degree polynomial that interpo-
lates it. We present a constructive proof of this fact by use of Lagrange Polynomials.
Definition 5.1 (Lagrange Polynomials). For a given set of n + 1 nodes Xj, the Lagrange
polynomials are the n + 1 polynomials £i defined by
(■i (Xj )
ifi^j
1 if i = j
Then we define the interpolating polynomial as
n
Pn{x) = ^2ydi(x).
i=0
If each Lagrange Polynomial is of degree at most n, then p n also has this property. The Lagrange
Polynomials can be characterized as follows:
ti(*)= n
(5.1)
63
64
CHAPTER 5. INTERPOLATION
By evaluating this product for each Xj , we see that this is indeed a characterization of the Lagrange
Polynomials. Moreover, each polynomial is clearly the product of n monomials, and thus has degree
no greater than n.
This gives the theorem
Theorem 5.2 (Interpolant Existence and Uniqueness). Let {xj}™ =0 be distinct nodes. Then
for any values at the nodes, {yi}™ =0 , there is exactly one polynomial, p(x) of degree no greater than
n such that p(xi) = yi for i = 0, 1, . . . , n.
Proof. The Lagrange Polynomial construction gives existence of such a polynomial p(x) of degree
no greater than n.
Suppose there were two such polynomials, call them p(x),q(x), each of degree no greater than
n, both interpolating the data. Let r(x) = p(x) — q(x). Note that r(x) can have degree no greater
than n, yet it has roots at xq, x\, . . . , x n . The only polynomial of degree < n that has n + 1 distinct
roots is the zero polynomial, i.e., = r(x) = p(x) — q(x). Thus p(x),q(x) are equal everywhere,
i.e., they are the same polynomial. □
Example Problem 5.3. Construct the polynomial interpolating the data
X
l
l
2
3
y
3
-10
2
by using Lagrange Polynomials.
Solution: We construct the Lagrange Polynomials:
e (x)
h(x)
Then the interpolating polynomial, in "Lagrange Form" is
p 2 (x) = -3(x - l -)(x - 3) - 8(x - l)(x - 3) + 2 -{x - l)(x - \)
H
5.1.2 Newton's Method
There is another way to prove existence. This method is also constructive, and leads to a different
algorithm for constructing the interpolant. One way to view this construction is to imagine how
one would update the Lagrangian form of an interpolant. That is, suppose some data were given,
and the interpolating polynomial calculated using Lagrange Polynomials; then a new point was
given (x n+ i, y n+ i), and an interpolant for the augmented set of data is to be found. Each Lagrange
Polynomial would have to be updated. This could take a lot of calculation (especially if n is large).
So the alternative method constructs the polynomials iteratively. Thus we create polynomials
Pk(x) such that Pk{xi) = y% for < i < k. This is simple for k = 0, we simply let
(x-i)(x-3)
(1-
-3)
(x —
l)(x
-3)
a-
-3)
(x —
l)(x
(3-l)(3-i)
~(x ~ ^)(x - 3)
|(* -!)(*-!)
Po(x) = y ,
5.1. POLYNOMIAL INTERPOLATION
65
2
-0.5 -
0.5 -
1.5 -
-1
-
1 -
0.5
1.5
2
2.5
3
3.5
Figure 5.1: The 3 Lagrange polynomials for the example are shown. Verify, graphically, that these
polynomials have the property ii(xj) = 5ij for the nodes 1, |,3.
the constant polynomial with value yo- Then assume we have a proper pk(x) and want to construct
Pk+i(x). The following construction works:
for some constant c. Note that the second term will be zero for any Xi for < i < k, so pk+i(x)
will interpolate the data at xq, x±, . . . , x^. To get the value of the constant we calculate c such that
Vk+i = Pk+i{xk+i) = Pk(xk+i) + c(xfc+i - x )(x k +i - xi) ■ ■ ■ - x k ).
This construction is known as Newton's Algorithm, and the resultant form is Newton's form of
the interpolant
Example Problem 5.4. Construct the polynomial interpolating the data
by using Newton's Algorithm.
Solution: We construct the polynomial iteratively:
p (x) = 3
Pi(x) = 3 + c(x — 1)
We want -10 = pi{\) = 3 + c(-|), and thus c = 26. Then
Pk+i{x) = Pk(x) + c(x - x )(x - xx) ■ ■ ■ (x - x k ),
x 1 A 3
y 3 -10 2
66
CHAPTER 5. INTERPOLATION
0.5 1 1.5 2 2.5 3 3.5 4
Figure 5.2: The interpolating polynomial for the example is shown.
We want 2 = p 2 (3) = 3 + 26(2) + c(2)(§), and thus c = =j£. Then we get
—53 1
p 2 (x) = 3 + 26(x - 1) + -— (x - l)(x - -).
5 2
H
Does Newton's Algorithm give a different polynomial? It is easy to show, by induction, that
the degree of p n (x) is no greater than n. Then, by Theorem 5.2, it must be the same polynomial
as the Lagrange interpolant.
The two methods give the same interpolant, we may wonder "Which should we use?" Newton's
Algorithm seems more flexible-it can deal with adding new data. There also happens to be a way
of storing Newton's form of the interpolant that makes the polynomial simple to evaluate (in the
sense of number of calculations required).
5.1.3 Newton's Nested Form
Recall the iterative construction of Newton's Form:
Pk+i(x) = Pk(x) + c k (x - x )(x -xi)---(x- x k ).
The previous iterate p k (x) was constructed similarly, so we can write:
Vk+i{x) = \Pk-i{x) + c fc _i(x - x )(x - Xi) ■ ■ ■ (x - x k _i)] + c k (x - x )(x -xi)---(x- x k ).
Continuing in this way we see that we can write
n
Pn(x) = y^ ck
k=0
0<j<k
where an empty product has value 1 by convention. This can be rewritten in a funny form, where
a monomial is factored out of each successive summand:
p n (x) = C + (X - X )[CI + {X- Xl)[c 2 + (x - X 2 ) [. . .]]]
5.1. POLYNOMIAL INTERPOLATION
67
Supposing for an instant that the constants c& were known, this provides a better way of
calculating p n (t) at arbitrary t. By "better" we mean requiring few multiplications and additions.
This nested calculation is performed iteratively:
V 1 = C n _i + {t- Xn^Vo
V2 = Cn-2 + (t - X n -2)Vl
V n = Co + (t- x )v n -i
This requires only n multiplications and 2n additions. Compare this with the number required
for using the Lagrange form: at least n 2 additions and multiplications.
5.1.4 Divided Differences
It turns out that the coefficients Ck for Newton's nested form can be calculated relatively easily by
using divided differences. We assume, for the remainder of this section, that we are considering
interpolating a function, that is, we have values of f(xi) at the nodes Xj.
Definition 5.5 (Divided Differences). For a given collection of nodes {xj}" =0 and values
{f( x i)}x=oi a ^ th order divided difference is a function of k + 1 (not necessarily distinct) nodes,
written as
/ [ x ioi x ii > • • • > x ik\
The divided differences are defined recursively as follows:
• The th order divided differences are simply defined:
f[Xi] = f(Xi).
• Higher order divided differences are the ratio of differences:
p r I / [ x ii i x i2 i • • • i x i}~\ ~ f \_ x io i x ii ) • • • )
/ [ x io ) x h > ■ ■ ■ > x ik\
We care about divided differences because coefficients for the Newton nested form are divided
differences:
Ck = f [x ,xi,.. .,x k \.
(5.3)
Because we are only interested in the Newton method coefficients, we will only consider divided
differences with successive nodes, i.e., those of the form / [xj, Xj+i, . . . , Xj+fc]. In this case the
higher order differences can more simply be written as
/ [Xj, Xj+l, . . . , Xj+k] —
f [ x j+l, x j+2, x j+k] ~ f [Xj,Xj + l, ... , X j+k _{\
x j+k x j
68
CHAPTER 5. INTERPOLATION
The graphical way to calculate these things is a "pyramid scheme" 1 , where you compute the
following table columnwise, from left to right:
X
f[]
/[,]
/[,,]
/[,,,]
x
f[xo]
f[xo,Xi]
X\
f[xi]
f[xo,xi
x 2 ]
f[xi,x 2 ]
f[x ,Xl,X 2 ,X3,]
X2
f[x2]
f[xi,x 2
X3]
f[x2,X 2 ]
X3
/N
Note that by definition, the first column just echoes the data: f[xi] = f{xi).
We drag out our example one last time:
Example Problem 5.6. Find the divided differences for the following data
X
1
1
2
3
Si?)
3
-10
2
Solution: We start by writing in the data:
X
/[]
/[,]
/[,,]
1
3
1
2
-10
3
2
Then we calculate:
f[xo,xi
(1/2) - 1
Adding these to the table, we have
2d, and j[xi,x 2 \
3 - (1/2)
24/5.
X
/[]
/[,]
1
3
26
1
2
-10
24
5
3
2
/[,,
Then we calculate
So we complete the table:
f[x ,xi,x 2 ]
24/5 - 26 _ -53
3-1 ~ ~~5~ '
X
/[]
/[,]
/[,,]
1
3
26
1
2
-10
-53
5
24
5
3
2
1 That's supposed to be a joke.
5.2. ERRORS IN POLYNOMIAL INTERPOLATION
69
You should verify that along the top line of this pyramid you can read off the coefficients for
Newton's form, as found in Example Problem 5.4. H
5.2 Errors in Polynomial Interpolation
We now consider two questions:
1. If we want to interpolate some function f(x) at n + 1 nodes over some closed interval, how
should we pick the nodes?
2. How accurate can we make a polynomial interpolant over a closed interval?
You may be surprised to find that the answer to the first question is not that we should make
the Xi equally spaced over the closed interval, as the following example illustrates.
Example 5.7 (Runge Function). Let
f(x) = {l + x 2 )-\
(known as the Runge Function), and let p n (x) interpolate / on n equally spaced nodes, including
the endpoints, on [—5,5]. Then
lim max \p n (x) — f(x)\ = oo.
n->oo x e[-5,5]
This behaviour is shown in Figure 5.3.
It turns out that a much better choice is related to the Chebyshev Polynomials ("of the first
kind"). If our closed interval is [—1,1], then we want to define our nodes as
Xi = cos
2i + l
2n + 2
7T
< i < n.
Literally interpreted, these Chebyshev Nodes are the projections of points uniformly spaced on
a semi circle; see Figure 5.4. By using the Chebyshev nodes, a good polynomial interpolant of the
Runge function of Example 5.7 can be found, as shown in Figure 5.5.
5.2.1 Interpolation Error Theorem
We did not just invent the Chebyshev nodes. The fact that they are "good" for interpolation
follows from the following theorem:
Theorem 5.8 (Interpolation Error Theorem). Let p be the polynomial of degree at most n
interpolating function / at the n+1 distinct nodes xq, xi, . . . , x n on [a, b]. Let /( n+1 ) be continuous.
Then for each x £ [a, b] there is some £ G [a, b] such that
1 "
/(*) " V{x) = 1 ——fW (0 H (x - Xl ) .
(n + 1)!
You should be thinking that the term on the right hand side resembles the error term in Taylor's
Theorem.
70
CHAPTER 5. INTERPOLATION
0.8
0.6
0.4
0.2
-0.2
-0.4
runge
3 nodes
7 nodes
10 nodes
-2 2
(a) 3,7,10 equally spaced nodes
(b) 15 equally spaced nodes
(c) 50 equally spaced nodes
Figure 5.3: The Runge function is poorly interpolated at n equally spaced nodes, as n gets large.
At first, the interpolations appear to improve with larger n, as in (a). In (b) it becomes apparent
that the interpolant is good in center of the interval, but bad at the edges. As shown in (c), this
gets worse as n gets larger.
5.2. ERRORS IN POLYNOMIAL INTERPOLATION
71
Figure 5.4: The Chebyshev nodes are the projections of nodes equally spaced on a semi circle.
(a) 3,7,10 Chebyshev nodes (b) 17 Chebyshev nodes
Figure 5.5: The Chebyshev nodes yield good polynomial interpolants of the Runge function. Com-
pare these figures to those of Figure 5.3. For more than about 25 nodes, the interpolant is indis-
tinguishable from the Runge function by eye.
Proof. First consider when x is one of the nodes xf, in this case both the LHS and RHS are zero.
So assume x is not a node. Make the following definitions
n
w{t) = n (t - ,
= f{x) -p(x)
w(x)
0(t) = f(t)-p(t)-cw(t).
Since x is not a node, w(x) is nonzero. Now note that 4>(xi) is zero for each node Xi, and that by
definition of c, that 4>(x) = for our x. That is 4>(t) has n + 2 roots.
Some other facts: f,p,w have n + 1 continuous derivatives, by assumption and definition; thus
cf) has this many continuous derivatives. Apply Rolle's Theorem to cj)(t) to find that (j)'(t) has n + 1
72
CHAPTER 5. INTERPOLATION
roots. Then apply Rolle's Theorem again to find 4>"(t) has n roots. In this way we see that (f) (n+l \t)
has a root, call it £. That is
o = 4> {n+1) (0 = f (n+1) (0-p {n+1) (0 -
But p(t) is a polynomial of degree < n, so j/ ra+1 ) is identically zero. And w(t) is a polynomial
of degree n + 1 in t, so its n + 1 th derivative is easily seen to be (n + 1)! Thus
= /(" +1 )(e)-c(n + l)!
c(n + l)! = /(" +1 )(e)
/(a:) -p(x)
w(x)
1
/(»+!)(£)
(n + 1)!'
which is what was to be proven.
Thus the error in a polynomial interpolation is given as
□
1 n
/(*) " V{x) = j—-^fW (£) H (x - x t ) .
We have no control over the function f(x) or its derivatives, and once the nodes and / are fixed, p
is determined; thus the only way we can make the error \f(x) — p(x)\ small is by judicious choice
of the nodes Xj.
The Chebyshev nodes on [—1, 1] have the remarkable property that
^=0
< 2"
for any t £ [—1, 1] . Moreover, it can be shown that for any choice of nodes Xi that
max
*e[-i,i]
II<«
> 2~ n .
Thus the Chebyshev nodes are considered the best for polynomial interpolation.
Merging this result with Theorem 5.8, the error for polynomial interpolants defined on Cheby-
shev nodes can be bounded as
I/O) ~p(x)\ <
max
2 n (n + l)Ue[-i,i]
f (n+1 Ho
The Chebyshev nodes can be rescaled and shifted for use on the general interval [a, (3\ . In this
case they take the form
(3 — a
Xi = cos
2i + l
2n + 2
7T
a + (3
0<i<n.
In this case, the rescaling of the nodes changes the bound on T\(t — %i) so the overall error bound
becomes
08 - «) n+1
\f(x)-p(x)\ <
max
2 2n+l ( n+ l)! 5e[Qj/3]
f(n + l) {0
for x £ [a, (3] .
5.2. ERRORS IN POLYNOMIAL INTERPOLATION
73
Example Problem 5.9. How many Chebyshev nodes are required to interpolate the function
f(x) = sin(x) + cos(x) to within 10~ 8 on the interval [0, 7r]?
Solution: We first find the derivatives of /
f'(x) = cos(x) — sin(x)
/ (x) = — sin(x) — cos(x)
f"'(x) = — cos(x) + sin(x)
As a crude approximation we can assert that
< Icosxl + I sin x I < 2.
Thus it suffices to take n large enough such that
:2 < 10- 8
2 2n+l ( n+ 1)!'
By trial and error we see that n = 10 suffices. H
5.2.2 Interpolation Error for Equally Spaced Nodes
Despite the proven superiority of Chebyshev Nodes, and the problems with the Runge Function,
equally spaced nodes are frequently used for interpolation, since they are easy to calculate 2 . We
now consider bounding
n
max TT |x — Xjl ,
where
Xi = a + hi = a + — — —i, i = 0, 1, . . . , n.
n
Start by picking an x. We can assume x is not one of the nodes, otherwise the product in
question is zero. Then x is between some Xj, Xj + \ We can show that
< — .
by simple calculus.
Now we claim that \x — Xj| < (j — i + 1) h for i < j, and |x — x«| < (i — j) h for j + 1 < i. Then
n h 2
Y[\x-Xi\< — [(j + 1)W] [(n - jW'i- 1 ] .
It can be shown that (j + l)!(n — j)\ < n\, and so we get an overall bound
h n+1 n\
X — X{\ <
i=0
Y\_ \X - Xi\ <
2 Never underestimate the power of laziness.
74
CHAPTER 5. INTERPOLATION
The interpolation theorem then gives us
h n+1
\f(x) —p(x)\ < —, — max
f(n + l) {0
where h = (b — a)/n.
The reason this result does not seem to apply to Runge's Function is that for Runge's
Function becomes unbounded as n — > oo.
Example Problem 5.10. How many equally spaced nodes are required to interpolate the function
f(x) = sin(x) + cos(a;) to within 10 -8 on the interval [0, 7r]?
Solution: As in Example Problem 5.9, we make the crude approximation
/<*>(*)
< Icosxl + I sin a; I < 2.
Thus we want to make n sufficiently large such that
h n+i
-2 < 10"°,
4(n+ 1)
where h = (tt — 0)/n. That is we want to find n large enough such that
7T ri+1
2n n+1 (n + 1)
By simply trying small numbers, we can see this is satisfied if n = 12.
5.2. ERRORS IN POLYNOMIAL INTERPOLATION
75
Exercises
(5.1) Find the Lagrange Polynomials for the nodes { — 1, 1}.
(5.2) Find the Lagrange Polynomials for the nodes { — 1, 1, 5}.
(5.3) Find the polynomial of degree no greater than 2 that interpolates
X
-l
1
5
y
3
3
-2
(5.4) Complete the divided differences table:
X
/[]
/[,]
/[,,]
-1
3
1
3
5
-2
Find the Newton form of the polynomial interpolant.
(5.5) Find the polynomial of degree no greater than 3 that interpolates
X
-l
1
5
-3
y
3
3
-2
4
(Hint: reuse the Newton form of the polynomial from the previous question.)
(5.6) Find the nested form of the polynomial interpolant of the data
X
l
3
4
6
y
-3
13
21
1
by completing the following divided differences table:
X
/[]
/[,]
/[,,]
/[,,,]
1
-3
3
13
4
21
6
1
(5.7) Find the polynomial of degree no greater than 3 that interpolates
X
l
3/2
2
y
3
2
37/8
8
(5.8) Complete the divided differences table:
X
/[]
/[,]
/[,,]
/[,,,]
-1
6
3
1
2
2
3
76
CHAPTER 5. INTERPOLATION
(Something a little odd should have happened in the last column.) Find the Newton form of
the polynomial interpolant. Of what degree is the polynomial interpolant?
(5.9) Let p(x) interpolate the function cos(x) at n equally spaced nodes on the interval [0, 2] . Bound
the error
max \p(x) — cos(x)|
0<x<2
as a function of n. How small is the error when n = 10? How small would the error be if
Chebyshev nodes were used instead? How about when n = 10?
(5.10) Let p(x) interpolate the function x~ 2 at n equally spaced nodes on the interval [0.5,1].
Bound the error
max \p(x) — x~ 2 \
0.5<a;<l 1 V ' 1
as a function of n. How small is the error when n = 10?
(5.11) How many Chebyshev nodes are required to interpolate the function ^ to within 10~ 6 on
the interval [1,2]?
(5.12) Write code to calculate the Newton form coefficients, by divided differences, for the nodes
Xi and values /(xj). Your m-file should have header line like:
function coefs = newtonCoef (xs, fxs)
where xs is the vector of n + 1 nodes, and f xs the vector of n + 1 values. Test your code on
the following input:
octave :1> xs = [1-12-23-34 -4] ;
octave :2> fxs = [1 1 2 3 5 8 13 21] ;
octave :3> newtonCoef (xs, fxs)
ans =
1.00000 -0.00000 0.33333 -0.08333 0.05000 0.00417 -0.00020 -0.00040
(a) What do you get when you try the following?
octave :4> xs = [1 4 5 9 14 23 37 60];
octave :5> fxs = [3 141592 6];
octave :6> newtonCoef (xs ,ys)
(b) Try the following:
octave :7> xs = [1 3 4 2 8 -2 14 23 15];
octave :8> fxs = xs.*xs + xs .+4;
octave :9> newtonCoef (xs, fxs)
(5.13) Write code to calculate a polynomial interpolant from its Newton form coefficients and the
node values. Your m-file should have header line like:
function y = calcNewton(t , coefs ,xs)
where coefs is a vector of the Newton coefficients, xs is a vector of the nodes Xj, and y is
the value of the interpolating polynomial at t. Check your code against the following values:
octave :1> xs = [13 4];
octave :2> coefs = [6 5 1];
octave :3> calcNewton (5,coefs,xs)
ans = 34
octave :4> calcNewton (-3 , coef s ,xs)
ans = 10
(a) What do you get when you try the following?
octave :5> xs= [3 14592687 0];
5.2. ERRORS IN POLYNOMIAL INTERPOLATION
77
octave :6> coefs = [1-12-23-34-45 -5];
octave : 7> calcNewton (0 . 5 , coefs , xs)
(b) Try the following
octave :8> xs= [3 14592687 0];
octave :9> coefs = [1-12-23-34-45 -5];
octave : 10> calcNewton(l , coef s ,xs)
CHAPTER 5. INTERPOLATION
Chapter 6
Spline Interpolation
Splines are used to approximate complex functions and shapes. A spline is a function consisting of
simple functions glued together. In this way a spline is different from a polynomial interpolation,
which consists of a single well defined function that approximates a given shape; splines are normally
piecewise polynomial.
6.1 First and Second Degree Splines
Splines make use of partitions, which are a way of cutting an interval into a number of subintervals.
Definition 6.1 (Partition). A partition of the interval [a, b] is an ordered sequence {tj}™ =0 such
that
a = t < h < ■ ■ ■ < t n -i <t n = b
The numbers ti are known as knots.
A spline of degree 1, also known as a linear spline, is a function which is linear on each subinterval
defined by a partition:
Definition 6.2 (Linear Splines). A function S is a spline of degree 1 on [a, b] if
1. The domain of S is [a, b\.
2. S is continuous on [a, b\.
3. There is a partition {ii}™ =0 of [a,b] such that on each [ti, , S is a linear polynomial.
A linear spline is defined entirely by its value at the knots. That is, given
t
to
h
y
yo
yi
y n
there is only one linear spline with these values at the knots and linear on each given subinterval.
For a spline with this data, the linear polynomial on each subinterval is defined as
Vi+i - yi
Si(x) = yi +
U+l ~ ti
(x - ti)
Note that if x € [U, U+i] , then x — ti > 0, but x — tj_i < 0. Thus if we wish to evaluate S(x), we
search for the largest i such that x — ti > 0, then evaluate Si(x).
Example 6.3. The linear spline for the following data
t
0.0
0.1
0.4
0.5
0.75
1.0
y
1.3
4.5
2.0
2.1
5.0
3
is shown in Figure 6.1.
79
80
CHAPTER 6. SPLINE INTERPOLATION
6
5 -
4 -
3 -
2 -
1 1 1 1 1 1 1 1
0.2 0.4 0.6 0.8 1
Figure 6.1: A linear spline. The spline is piecewise linear, and is linear between each knot.
6.1.1 First Degree Spline Accuracy
As with polynomial functions, splines are used to interpolate tabulated data as well as functions. In
the latter case, if the spline is being used to interpolate the function /, say, then this is equivalent
to interpolating the data
t
to
h
tn
V
f(to)
f(tn)
A function and its linear spline interpolant are shown in Figure 6.2. The spline interpolant in
that figure is fairly close to the function over some of the interval in question, but it also deviates
greatly from the function at other points of the interval. We are interested in finding bounds on
the possible error between a function and its spline interpolant.
To find the error bound, we will consider the error on a single interval of the partition, and use
a little calculus. 1 Suppose p(t) is the linear polynomial interpolating f(t) at the endpoints of the
subinterval then for t G [ij,tj + i] ,
\f(t) - p(t)\ < max{|/(i) - f(U)\ , \f(t) - f(t i+ i)\} .
That is, \f(t) — p(t)\ is no larger than the "maximum variation" of f(t) on this interval.
In particular, if f'(t) exists and is bounded by M\ on [if, tj+i], then
\f(t)-p(t)\<^(t i+1 -U).
Similarly, if f"(t) exists and is bounded by M2 on then
\f(t)-p(t)\<^(t l+1 - ti ) 2 .
Over a given partition, these become, respectively
Although we could directly claim Theorem 5.8, it is a bit of overkill.
6.1. FIRST AND SECOND DEGREE SPLINES
81
Figure 6.2: The function f(t) = 4.1 + sin (1/(0. 08i + 0.5)) is shown, along with the linear spline
interpolant, for the knots 0,0.1,0.4,0.6,0.75,0.9,1.0 For some values of t, the function and its
interpolant are very close, while they vary greatly in the middle of the interval.
-p(t)\ <
Mi
~Y
max (ti+i -U) ,
0<i<n
\f(t)
-v{t)\ <
M 2
8
max (t i+1 - ti) 2 .
0<i<n
(6.1)
If equally spaced nodes are used, these bounds guarantee that spline interpolants become better
as the number of nodes is increased. This contrasts with polynomial interpolants, which may get
worse as the number of nodes is increased, cf. Example 5.7.
6.1.2 Second Degree Splines
Piecewise quadratic splines, or splines of degree 2, are defined similarly:
Definition 6.4 (Quadratic Splines). A function Q is a quadratic spline on [a,b] if
1. The domain of Q is [a,b].
2. Q is continuous on [a, b].
3. Q' is continuous on (a,b).
4. There is a partition {ti}" =0 of [a, b] such that on Q is a polynomial of degree at most
2.
Example 6.5. The following is a quadratic spine:
Q(x) =
x
—x x < 0,
2 x < x < 2,
-x 2 + 7x-8 2<x.
Unlike linear splines, quadratic splines are not defined entirely by their values at the knots. We
consider why that is. The spline Q{x) is defined by its piecewise polynomials,
Qi(x) = diX 2 + biX + Cj.
82
CHAPTER 6. SPLINE INTERPOLATION
Thus there are 3n parameters to define Q(x).
For each of the n subintervals, the data
t
to
h
V
yo
yi
Vn
give two equations regarding Qi(x), namely that Qi(U) must equal yi and must equal j/i+i-
This is 2n equations. The condition on continuity of Q' gives a single equation for each of the n — 1
internal nodes. This totals 3n — 1 equations, but 3n unknowns. This system is under determined.
Thus some additional user-chosen condition is required to determine the quadratic spline. One
might choose, for example, Q'(a) = 0, or Q"(a) = 0, or some other condition.
6.1.3 Computing Second Degree Splines
Suppose the data
t
to
h
tn
y
yo
yi
y n
are given. Let Zi = Q'^U), and suppose that the additional condition to define the quadratic spline
is given by specifying zq. We want to be able to compute the form of Qi(x).
Because Qi(U) = yi, Q'^U) = Zi, Q'^U+i) = we see that we can define
Qi( x ) = wrr- — tt ( x ~ + z i(x~ U) + yi.
* {H+1 — H)
Use this at tj+i :
y i+ i = Qi(t i+ i) = — ^— (t i+ i - ti) 2 + Zi (t i+ i - U) + y i}
* {H+l — ti)
yi+i -yi = Zt+1 2 (U+i - ti) + Zi (t i+ i - u) ,
Vi+i ~Vi = 2 ~ tl > ■
Thus we can determine, from the data alone, Zi + \ from Zf.
_ n yi+i - yi
Zi+l — Zi.
H+l — ti
6.2 (Natural) Cubic Splines
If you recall the definition of the linear and quadratic splines, probably you can guess the definition
of the spline of degree k:
Definition 6.6 (Splines of Degree k). A function S is a spline of degree k on [a, b] if
1. The domain of S is [a, b\.
2. S, S', S", S^- 1 ) are continuous on (a, b).
3. There is a partition {tj}™ =0 of [a, b] such that on [ti, U+i], S is a polynomial of degree < k.
You would also expect that a spline of degree k has k — 1 "degrees of freedom," as we show
here. If the partition has n + 1 knots, the spline of degree k is defined by n(k + 1) parameters. The
given data
t
to
h
tn
y
yo
yi
y n
6.2. (NATURAL) CUBIC SPLINES
83
provide 2n equations. The continuity of S', S" , . . . , S^^ 1 ^ at the n — 1 internal knots gives (k —
l)(n — 1) equations. This is a total of n(k + 1) — (k — 1) equations. Thus we have k — 1 more
unknowns than equations. Thus, barring some singularity, we can (and must) add k — 1 constraints
to uniquely define the spline. These are the degrees of freedom.
Often k is chosen as 3. This yields cubic splines. We must add 2 extra constraints to define the
spline. The usual choice is to make
S"(to) = S"(t n ) = 0.
This yields the natural cubic spline.
6.2.1 Why Natural Cubic Splines?
It turns out that natural cubic splines are a good choice in the sense that they are the "inter-
polant of minimal H 2 seminorm." The corollary following this theorem states this in more easily
understandable terms:
Theorem 6.7. Suppose / has two continuous derivatives, and S is the natural cubic spline inter-
polating / at knots a = to < t\ < . . . < t n = b. Then
f [S"{x)] 2 dx < f [f"{x)] 2 dx
J a J a
Proof. We let g{x) = f(x) — S(x). Then g(x) is zero on the (n + 1) knots t{. Derivatives are linear,
meaning that
f"(x) = S"(x)+g"(x).
Then b b b b
f [f"(x)] 2 dx= f [S"(x)] 2 dx+ f [g"(x)] 2 dx+ f 2S"(x)g"(x)dx.
J a J a J a J a
We show that the last integral is zero. Integrating by parts we get
f S"(x)g"(x)dx = S"g'
J a
because S"(a) = S"(b) = 0. Then notice that S is a polynomial of degree < 3 on each interval, thus
S"'(x) is a piecewise constant function, taking value Cj on each interval [ij,tj + i]. Thus
n— 1 „h n—1
r-b b pb rb
S"'g'dx = - / S'"g'dx,
/ s'" g ' dx = Y, c l9 ' dx = ag\ = o,
i=0 " u i=0
with the last equality holding because g(x) is zero at the knots. □
Corollary 6.8. The natural cubic spline is best twice-continuously differentiable interpolant for a
twice-continuously differentiable function, under the measure given by the theorem.
Proof. Let / be twice-continuously differentiable, and let S be the natural cubic spline interpolating
f(x) at some given nodes {ij}" =0 . Let R(x) be some twice-continuously differentiable function
which also interpolates f(x) at these nodes. Then S(x) interpolates R(x) at these nodes. Apply
the theorem to get
f [S"(x)} 2 dx< f [R"{x)] 2 dx
J a J a
□
84
CHAPTER 6. SPLINE INTERPOLATION
6.2.2 Computing Cubic Splines
First we present an example of computing the natural cubic spline by hand:
Example Problem 6.9. Construct the natural cubic spline for the following data:
t
-l
2
y
3
-1
3
Solution: The natural cubic spline is defined by eight parameters:
S(x) =
We interpolate to find that d = h = — 1 and
ax 3 + bx 2 + cx + d x G [-1, 0]
ex 3 + fx 2 + gx + h x G [0, 2]
We take the derivative of S:
S'(x) =
—a + 6 — c — 1 = 3
8e + 4/ + 2g - 1 = 3
3ax 2 + 2bx + c x G [-1,0]
3ex 2 + 2fx + g x G [0, 2]
Continuity at the middle node gives c = g. Now take the second derivative of 5:
„ f 6ax + 2b x G [-1,0]
lX j ~ \ 6ex + 2/ x G [0, 2]
Continuity at the middle node gives b = f. The natural cubic spline condition gives —6a + 2b =
and 12e + 2f = 0. Solving this by "divide and conquer" gives
o/ n J ar 3 + 3x 2 - 2x - 1 xG[-l,0]
(Xj ~ \ -\x 3 + 3x 2 - 2x - 1 x G [0, 2]
H
Finding the constants for the previous example was fairly tedious. And this is for the case of
only three nodes. We would like a method easier than setting up the \n equations and unknowns,
something akin to the description in Subsection 6.1.3. The method is rather tedious, so we leave it
to the exercises.
6.3 B Splines
The B splines form a basis for spline functions, whence the name. We presuppose the existence of
an infinite number of knots:
. . . < t 2 < h < t < h < t 2 < ■ ■ ■ ,
with lim^-oo tk = —oo and lim^oo tk = oo.
The B splines of degree are defined as single "blocks" :
S?(x)
1 ti < X < t i+ l
otherwise
6.3. B SPLINES
85
The zero degree B splines are continuous from the right, are nonzero only on one subinterval
[U,ti + i), sum to 1 everywhere.
We justify the description of B splines as basis splines: If S is a spline of degree on the given
knots and is continuous from the right then
That is, the basis splines work in the same way that Lagrange Polynomials worked for polynomial
interpolation.
The B splines of degree k are defined recursively:
Some B splines are shown in Figure 6.3.
The B splines quickly become unwieldy. We focus on the case k = 1. The B spline Bj(x) is
• Piecewise linear.
• Continuous.
• Nonzero only on (^,£3+2).
These B splines are sometimes called hat functions. Imagine wearing a hat shaped like this! What-
ever.
The nice thing about the hat functions is they allow us to use analogy. Harken back to polyno-
mial interpolation and the Lagrange Functions. The hat functions play a similar role because
S(x) = Y,S(x l )B?(x).
• 1 at ij+i.
Then if we want to interpolate the following data with splines of degree 1:
t to t\
t
n
v yo vi
We can immediately set
n
s(x) = Y,y l Bl 1 {x).
i=0
86
CHAPTER 6. SPLINE INTERPOLATION
Figure 6.3: Some B splines for the knots to = 1, ti = 2.5, £2 = 3.5,^3 = 4 are shown. In (a), one of
the degree B splines; in (b), a degree 1 B spline; in (c), two of the degree 1 B splines and a degree
2 B spline are shown. The two hat functions "merge" together to give the quadratic B spline.
6.3. B SPLINES
87
Exercises
S(x) =
S(x)
on [0,4]?
Why or why
~3x + 2
: < x < 1
-2x + 4
: 1 < x < 4
on [0,2]?
Why or why
fx + 3
: < x < 1
I 3
: 1 < x < 2
on [0,4]?
Why or why
fx 2 + 3
: < x < 3
1 5x — 6
: 3 < x < 4
(6.4) Find constants, a, [3 such that the following is a linear spline on [0,5].
S{x)
Ax -2
ax + (5
-2x + 10
< x < 1
1 < x < 3
3 < x < 5
(6.5) Is the following function a quadratic spline on [0, 4]? Why or why not?
, fx 2 + 3 : < x < 3
Q(x) = <
[5x-6 : 3 < x < 4
(6.6) Is the following function a quadratic spline on [0, 2]? Why or why not?
'x 2 + 3x + 2 :0<x<l
Q(x)
2x 2 + x + 3 :l<x<2
(6.7) Find constants, a, (3, 7 such that the following is a quadratic spline on [0,5].
Q(x) = <
(6.8) Find the quadratic spline that interpolates the following data:
'ix 2 + 2x + §
ax 2 + 0x + 7
3x 2 - 7x + 12
< x < 1
1 < x < 3
3 < x < 5
t
1
4
y
1
-2
1
To resolve the single degree of freedom, assume that Q'(0) = —Q'(4). Assume your solution
takes the form
|ai(x-l) 2 + /3i(x-l)-2 :0<x<l
|a 2 (x-l) 2 + /3 2 (x-l)-2 :l<x<4
Find the constants a±, Pi, a 2 , /?2-
Q(x)
88
CHAPTER 6. SPLINE INTERPOLATION
(6.9) Find the natural cubic spline that interpolates the data
X
1
3
y
4
2
7
It may help to assume your answer has the form
S(X):
Ax 3 + Bx 2 + Cx + 4
D(x - l) 3 + E(x - l) 2 + F(x-l) + 2
< x < 1
1< x < 3
Chapter 7
Approximating Derivatives
7.1 Finite Differences
Suppose we have some blackbox function f(x) and we wish to calculate f'{x) at some given x. Not
surprisingly, we start with Taylor's theorem:
f( x + h ) = f(x) + f'(x)h+^-.
Rearranging we get
_ f( x + h)-f(x) _ f"(Qh
nX) ~ h 2 '
Remember that £ is between x and x + h, but its exact value is not known. If we wish to calculate
f'(x), we cannot evaluate /"(£), so we approximate the derivative by dropping the last term. That
is, we calculate [f(x + h) — f(x)] /h as an approximation 1 to f'(x). In so doing, we have dropped
the last term. If there is a finite bound on f"(z) on the interval in question then the dropped term
is bounded by a constant times h. That is,
fix) = + + O (h)
(7.1)
The error that we incur when we approximate f'(x) by calculating [f(x + h) — f(x)] /his called
truncation error. It has nothing to do with the kind of error that you get when you do calculations
with a computer with limited precision; even if you worked in infinite precision, you would still
have truncation error.
The truncation error can be made small by making h small. However, as h gets smaller, precision
will be lost in equation 7.1 due to subtractive cancellation. The error in calculation for small h
is called roundoff error. Generally the roundoff error will increase as h decreases. Thus there is
a nonzero h for which the sum of these two errors is minimized. See Example Problem 7.5 for an
example of this.
The truncation error for this approximation is O (h). We may want a more precise approxima-
tion. By now, you should know that any calculation starts with Taylor's Theorem:
f(x + h) = f( x ) + f'( x )h + ^lh 2 + ^^h 3
f(x-h) = f(x)-f'(x)h + ^lh 2 -^Mh"
lr This approximation for f'(x) should remind you of the definition of f'(x) as a limit.
89
90
CHAPTER 7. APPROXIMATING DERIVATIVES
By subtracting these two lines, we get
f{x + h)- f(x - h)
Thus
2f'(x)h
f'(x)
2f\x)h+ f " { ^ )+ J" {i2) h Z .
f(x + h)- f{x - h)
fjx + h)- fjx - h)
2h
r«i)+r(6) fc s
3!
h 2
f ,n (P \4-f'"(P \
If f (x) is continuous, then there is some f between £ 1; £ 2 such that /"'(£) = 1 ^ 2> . (This
is the MVT at work.) Assuming some uniform bound on /'"(■), we get
/ , W = /( ' + lk) - /(x -* ) +o(tf)
(7.2)
In some situations it may be necessary to use evaluations of a function at "odd" places to
approximate a derivative. These are usually straightforward to derive, involving the use of Taylor's
Theorem. The following examples illustrate:
Example Problem 7.1. Use evaluations of / at x + h and x + 2h to approximate fix), assuming
fix) is an analytic function, i.e., one with infintely many derivatives.
Solution: First use Taylor's Theorem to expand f{x + h) and f{x + 2h), then subtract to get some
factor of fix):
fix + 2h)
fix + h)
fix) + 2hfix) + yg-fix) + ^f'ix) + ^f 4 \x
+
fix) + hfix) + % f\x) + %f"ix) + %f 4 \x) + ...
fix + 2h)-fix + h) = hfix) + ^fix)+ 7 -^f'ix) + 1 -^f 4 \x) + .
ifix + 2h)-fix + h))/h = f(x) + fr(x) + ^r(x) + ifi/( 4 )(x) + ...
Thus ifix + 2h) - fix + h)) /h = fix) + O ih)
Example Problem 7.2. Show that
4fix + h)-fix + 2h)-3fix)
2h
= fix) + 0(h 2 ),
for / with sufficient number of derivatives
Solution: In this case we do not have to find the approximation scheme, it is given to us. We only
have to expand the appropriate terms with Taylor's Theorem. As before:
fix + h) = /(x) + TO + f/"(x) + S/'"(x) + ...
Afix + h) = Afix) + Ahfix) + ^f"ix) + ^f"ix) + ...
fix + 2h) = fix) + 2hfix) + ^fix) + ^nx) + ...
4/(s + h)- fix + 2h)
Afix + h)- fix + 2h) - 3/(x)
i4fix + h)-fix + 2h)-3fix))/2h
3/(x) + 2hfix) + Ofix) + =^f"ix) +
2hf'[x) + + . .
f'(x) + =^-f'"(x) + ...
7.2. RICHARDSON EXTRAPOLATION
91
7.1.1 Approximating the Second Derivative
Suppose we want to approximate the second derivative of some blackbox function f(x). Again,
start with Taylor's Theorem:
f(x + h) = f(x) + f\x)h+^h 2 + ^h 3 + ^h 4 + ...
f(x-h) = f(x)-f\x)h + l^hi-^h* + i^h*-...
Now add the two series to get
.2^ , o/ (4) W,4 , / (6) (^),6
f(x + h) + f(x-h) = 2f(x) + h z f"{x) + 2 J -^h' i + 2 J —^h!> + ...
Then let
f( x + h)-2f(x) + f(x-h) t „. . , J^(x) u2 , jW(x) u4
m = - — ^2 " = f"{x) + r-^h 2 + 2 J -^h A + ...,
oo
2fc
fc=i
Thus we can use Richardson Extrapolation on ip(h) to get higher order approximations.
This derivation also gives us the centered difference approximation to the second derivative:
nx)= f(* + h)-2f(x) + f(x- h ) +o{h2y
(7.3)
7.2 Richardson Extrapolation
The centered difference approximation gives a truncation error of O (/i 2 ) , which is better than
O (h) . Can we do better? Let's define
m = ^[f(x+h)-f(x-h)}.
Had we expanded the Taylor's Series for f(x + h), f(x — h) to more terms we would have seen
that
0(h) = f{x) + a 2 h 2 + a 4 h 4 + a 6 h 6 + a s h s + ...
The constants ctj are a function of f^ +1 \x) only. (In fact, they should take the value of ■)
What happens if we now calculate 4>(h/2)7
<f>(h/2) = f'( x ) + ^a 2 h 2 + ^a 4 h 4 + ^a 6 h e + ^a s h 8 + ...
But we can combine this with (f>{h) to get better accuracy. We have to be a little tricky, but we
can get the O (/i 2 ) terms to cancel by taking the right multiples of these two approximations:
4,(h) - = -3f(x) + 1 4h* + -a fj k l> + -a a h* + ...
3 v 1 4 16 64
92
CHAPTER 7. APPROXIMATING DERIVATIVES
This approximation has a truncation error of O (/i 4 ) .
This technique of getting better approximations is known as the Richardson Extrapolation, and
can be repeatedly applied. We will also use this technique later to get better quadrature rules-that
is, ways of approximating the definite integral of a function.
7.2.1 Abstracting Richardson's Method
We now discuss Richardson's Method in a more abstract framework. Suppose you want to calculate
some quantity L, and have found, through theory, some approximation:
<j>{h) = L + a 2kh
2k
k=l
Let
Now define
D(n,0)
£(n, m)
4 m £(n, m - 1) - £(n - 1, m - 1)
4m _ J
(7.4)
We will be interested in calculating D(n,n) for some n. We claim that
D{n,n) = L + o(h 2(n+ Vy
First we examine the recurrence for £(n, m). As in divided differences, we use a pyramid table:
0(0,0)
0(1,0) 0(1,1)
£(2,0) £(2,1) £(2,2)
£(n,0) £(n, 1) £(ra,2)
£(n, n)
By definition we know how to calculate the first column of this table; every other entry in the table
depends on two other entries, one directly to the left, and the other to the left and up one space.
Thus to calculate £(n,n) we have to compute this whole lower triangular array.
We want to show that £(n, n) = L + O (/i 2 ( n+1 )) , that is £(n, n) is a O (/i 2(n+1 )) approximation
to L. The following theorem gives this result:
Theorem 7.3 (Richardson Extrapolation). There are constants ak, m such that
00 / h\ 2k
D(n,m) = L+ ^ a k ^ m ( — J (0 < m < n) .
k=m+l ^ '
The proof is by an easy, but tedious, induction. We skip the proof.
7.2. RICHARDSON EXTRAPOLATION
93
7.2.2 Using Richardson Extrapolation
We now try out the technique on an example or two.
Example Problem 7.4. Approximate the derivative of f(x) = logx at x = 1.
Solution: The real answer is /'(l) = 1/1 = 1, but our computer doesn't know that. Define
*W-^[«i + *)-/<i-*)l = ^.
Let's use h = 0.1. We now try to find D(2, 2), which is supposed to be a O (/i 6 ) approximation to
/'(I) = 1:
n\m
1
2
« 1.003353477
1
log ^l « 1.000834586
« 0.999994954
2
g l « 75 « 1.000208411
« 0.999999686
» 1.000000002
This shows that the Richardson method is pretty good. However, notice that for this simple
example, we have, already, that 0(0.00001) w 0.999999999. H
Example Problem 7.5. Consider the ugly function:
f{x) = arctan(x).
Attempt to find f'{y/2). Recall that fix) = j^i, so the value that we are seeking is |.
Solution: Let's use h = 0.01. We now try to find D(2, 2), which is supposed to be a O (/i 6 )
approximation to |:
n\m
1
2
0.333339506181068
1
0.333334876543723
0.333333333331274
2
0.33333371913582
0.333333333333186
0.333333333333313
Note that we have some motivation to use Richardson's method in this case: If we let
m = ^[f(V2+h)-f(V2-h)
then making h small gives a good approximation to /' (V2~) until subtractive cancelling takes over.
The following table illustrates this:
94
CHAPTER 7. APPROXIMATING DERIVATIVES
h
<Kh)
1.0
0.39269908169872408
0.1
0.33395069677431943
0.01
0.33333950618106845
0.001
0.33333339506169679
0.0001
0.33333333395058062
1 x 10"
-5
0.33333333334106813
1 x 10"
-6
0.33333333332441484
1 x 10"
-7
0.33333333315788138
1 x 10"
-8
0.33333332760676626
1 x 10"
-9
0.33333336091345694
1 x 10"
10
0.333333360913457
1 x lO -
11
0.333333360913457
1 x 10-
12
0.33339997429493451
1 x 10-
13
0.33306690738754696
1 x 10-
14
0.33306690738754696
1 x 10-
15
0.33306690738754691
1 x 10-
16
The data are illustrated in Figure 7.1. Notice that <j)(h) gives at most 10 decimal places of
accuracy, then begins to deteriorate; Note however, we get 13 decimal places from D(2,2). H
1e-12 1 ■ 1 ■ ' ■ ' ■ ' ■ ' ■ ' ■ ' ■ 1
1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1
Figure 7.1: The total error for the centered difference approximation to f'(y/2) is shown versus
h. The total error is the sum of a truncation term which decreases as h decreases, and a roundoff
term which increases. The optimal h value is around 1 x 10 -5 . Note that Richardson's D(2, 2)
approximation with h = 0.01 gives much better results than this optimal h.
7.2. RICHARDSON EXTRAPOLATION
95
Exercises
(7.1) Derive the approximation
/'(*)
4/ (x + h)- 3/(x) - fjx - 2h)
6h
(7.2)
using Taylor's Theorem.
(a) Assuming that f(x) has bounded derivatives, give the accuracy of the above approxi-
mation. Your answer should be something like O (h 7 ).
(b) Let f(x) = x 3 . Approximate /'(0) with this approximation, using h = \.
Let f(x) be an analytic function, i.e., one which is infinitely differentiable. Let ip(h) be the
centered difference approximation to the first derivative:
f( x + h)-f(x-h)
2h
(a) Show that 4>(h) = f'(x) + f f"'(x) + ^f^(x) + ^f^(x) + . .
(b) Show that
8 (« W -»/2|) =rw+0(jI)i
(7.3) Derive the approximation
h 2
4/(x + 3h) + 5f(x) - 9fjx - 2h)
30h
using Taylor's Theorem.
(a) What order approximation is this? (Assume f(x) has bounded derivatives of arbitrary
order.)
(b) Use this formula to approximate /'(0), where f(x) = x 4 , and h = 0.1
(7.4) Suppose you want to know quantity Q, and can approximate it with some formula, say (f)(h),
which depends on parameter h, and such that <p(h) = Q + a±h + a2h 2 + a^h 3 + a^h 4 + . . .
Find some linear combination of (f>(h) and <j>{—h) which is a O (/i 2 ) approximation to Q.
(7.5) Assuming that 4>(h) = Q + a2h 2 + a^h 4 + aeh 6 . . ., find some combination of 4>(h), 4>(h/3)
which is a O (/i 4 ) approximation to Q.
(7.6) Let A be some number in (0, 1). Assuming that <j)(h) = Q + a^l? + a^h 4 + agh 6 . . ., find some
combination of (f)(h) , 4>(Xh) which is a O (h 4 ) approximation to Q. To make the constant
associated with the h 4 term small in magnitude, what should you do with A? Is this practical?
Note that the method of Richardson Extrapolation that we considered used the value A = 1/2.
(7.7) Assuming that <p{h) = Q + a2h 2 + a^h 4 + a^h 6 . . ., find some combination of (p(h),(j>(h/4)
which is a O (h 4 ) approximation to Q.
(7.8) Suppose you have some great computational approximation to the quantity Q such that
ij)(h) = Q + a%h 3 + a^h 6 + agh 9 . . . Can you find some combination of ip(h),ip(h/2) which is
a O (/i 6 ) approximation to Ql
(7.9) Complete the following Richardson's Extrapolation Table, assuming the first column consists
of values D(n, 0) for n = 0, 1, 2:
n\m
1
2
2
1
1.5
?
2
1.25
?
?
(See equation 7.4 if you've forgotten the definitions.)
CHAPTER 7. APPROXIMATING DERIVATIVES
10) Write code to complete a Richardson's Method table, given the first column. Your m-file
should have header line like:
function Dnn = richardsons (colO)
where Dnn is the value at the lower left corner of the table, D(n, n) while colO is the column
of n + 1 values D(i, 0), for i = 0, 1, . . . , n. Test your code on the following input:
octave :1> colO = [1 0.5 0.25 0.125 0.0625 0.03125];
octave :2> richardsons (colO)
ans = 0.019042
(a) What do you get when you try the following?
octave :5> colO = [1.5 0.5 1.5 0.5 1.5 0.5 1.5];
octave :6> richardsons (colO)
(b) What do you get when you try the following?
octave :7> colO = [0.9 0.99 0.999 0.9999 0.99999];
octave :8> richardsons (colO)
Chapter 8
Integrals and Quadrature
8.1 The Definite Integral
Often enough the numerical analyst is presented with the challenge of finding the definite integral
of some function:
In your golden years of Calculus, you learned the Fundamental Theorem of Calculus, which claims
that if /(x) is continuous, and F(x) is an antiderivative of f(x), then
What you might not have been told in Calculus is there are some functions for which a closed
form antiderivative does not exist or at least is not known to humankind. Nevertheless, you may
find yourself in a situation where you have to evaluate an integral for just such an integrand. An
approximation will have to do.
8.1.1 Upper and Lower Sums
We will review the definition of the Riemann integral of a function. A partition of an interval [a, b]
is a finite, ordered collection of nodes xf.
Given such a partition, P, define the upper and lower bounds on each subinterval [
follows:
a = xq < x\ < X2 < ■ ■ ■ < x.
n
= b.
mi = inf {/(x) | Xi < x < x i+1 }
Mi = sup {/(x) | Xi < x < x i+1 }
Then for this function / and partition P, define the upper and lower sums:
n-1
i=0
n-1
U(f,P)
^2 Mi (x i+ i - x^
i=0
97
CHAPTER 8. INTEGRALS AND QUADRATURE
9N
We can interpret the upper and lower sums graphically as the sums of areas of rectangles defined
by the function / and the partition P, as in Figure 8.1.
(a) The Lower Sum
(b) The Upper Sum
Figure 8.1: The (a) lower, and (b) upper sums of a function on a given interval are shown. These
approximations to the integral are the sums of areas of rectangles. Note that the lower sums are
an underestimate, and the upper sums an overestimate of the integral.
Notice a few things about the upper, lower sums:
(ii) If we switch to a "better" partition (i.e., a finer one), we expect that L(/, •) increases and
U(f,-) decreases.
The notion of integrability familiar from Calculus class (that is Riemann Integrability) is defined
in terms of the upper and lower sums.
Definition 8.1. A function / is Riemann Integrable over interval [a, b] if
where the supremum and infimum are over all partitions of the interval [a, b\. Moreover, in case
f(x) is integrable, we define the integral
You may recall the following
Theorem 8.2. Every continuous function on a closed bounded interval of the real line is Riemann
Integrable (on that interval).
Continuity is sufficient, but not necessary.
Example 8.3. Consider the Heaviside function:
This function is not continuous on any interval containing 0, but is Riemann Integrable on every
closed bounded interval.
(i) L(f,P)<U(f,P).
supL(/,P)=inf£/(/,P)
p P
P
8. 1 . THE DEFINITE INTEGRAL
99
Example 8.4. Consider the Dirichlet function:
. . _ f x rational
} 1 i irrational
For any partition P of any interval [a, b], we have L(f, P) = 0, while U (/, P) = 1, so
supL(/,P) = 0^1 = inf£/(/,P),
p p
so this function is not Riemann Integrable.
8.1.2 Approximating the Integral
The definition of the integral gives a simple method of approximating an integral f(x) dx. The
method cuts the interval [a, b] into a partition of n equal subintervals X{ = a+^p, for i = 0, 1, . . . , n.
The algorithm then has to somehow find the supremum and infimum of f{x) on each interval
[xj,Xj + i]. The integral is then approximated by the mean of the lower and upper sums:
cb
J f{x)6xKt\{L{f,P) + U{f,P))
Because the value of the integral is between L(f,P) and U(f,P), this approximation has error at
most
\ (U(f, P)~L(f,P)).
Note that in general, or for a black box function, it is usually not feasible to find the suprema
and infima of f(x) on the subintervals, and thus the lower and upper sums cannot be calculated.
However, if some information is known about the function, it becomes easier:
Example 8.5. Consider for example, using this method on some function f(x) which is monotone
increasing, that is x < y implies f(x) < f(y). In this case, the infimum of f(x) on each interval
occurs at the leftmost endpoint, while the supremum occurs at the right hand endpoint. Thus for
this partition, P, we have
n— 1 I n— 1
£(/> P ) = ^2 m i \ x k+l ~ Xk\ = ^2 f( X k)
k=0 k=0
n—l |, | n— 1 i, I n
u(f, p) = J2^ - x k\ = E f( x ^ = E
k=0 k=0 k=l
Then the error of the approximation is
8.1.3 Simple and Composite Rules
For the remainder of this chapter we will study "simple" quadrature rules, i.e., rules which approx-
imate the integral of a function, f(x) over an interval [a, b] by means of a number of evaluations of
/ at points in this interval. The error of a simple quadrature rule usually depends on the function
100
CHAPTER 8. INTEGRALS AND QUADRATURE
f, and the width of the interval [a, b] to some power which is determined by the rule. That is we
usually think of a simple rule as being applied to a small interval.
To use a simple rule on a larger interval, we usually cast it into a "composite" rule. Thus the
trapezoidal rule, which we will study next becomes the composite trapezoidal rule. The means of
extending a simple rule to a composite rule is straightforward: Partition the given interval into
subintervals, apply the simple rule to each subinterval, and sum the results. Thus, for example if
the interval in question is [a, 0\, and the partition is a = xq < x\ < X2 < ■ ■ ■ < x n = j3, we have
n-1
composite rule on [a, (3] = simple rule applied to
i=0
8.2 Trapezoidal Rule
Suppose we are trying to approximate the integral
f{x) dx,
for some unpleasant or black box function f(x).
The trapezoidal rule approximates the integral
/ f{x) dx
J a
by the (signed) area of the trapezoid through the points (a, /(a)) , (6, f(b)) , and with one side the
segment from a to b. See Figure 8.2.
Figure 8.2: The trapezoidal rule for approximating the integral of a function over [a, b] is shown.
By old school math, we can find this signed area easily. This gives the (simple) trapezoidal rule:
8.2. TRAPEZOIDAL RULE
101
The composite trapezoidal rule can be written in a simplified form, one which you saw in your
calculus class, if the interval in question is partitioned into equal width subintervals. That is if you
let [a, P] be partitioned by
a = xq < x\ < X2 < x 3 < . . . < x n = (3,
with Xi = a + ih, where h = ((3 — a)/n, then the composite trapezoidal rule is
/ f{x) dx = £ / f(?) ~ ? E (^+1 " x *) \f( x *) + • t 8 - 1 )
1/0 i=0 ^ i=0
Since each subinterval has equal width, Xi + \ — x% = h, and we have
(8.2)
In your calculus class, you saw this in the less comprehensible form:
["fix)
J a
dx « h
f(a) + f(b)
n-l
i=i
Note that the composite trapezoidal rule for equal subintervals is the same as the approximation
we found for increasing functions in Example 8.5.
Example Problem 8.6. Approximate the integral
f—
Jo 1+:
dx
by the composite trapezoidal rule with a partition of equally spaced points, for n = 2.
Solution: We have h = = 1, and /(xo) = l,/(xi) = ^,/(x2) = \. Then the composite
trapezoidal rule gives the value
I[ / (x ) + /(x 1 ) + /(x 1 ) + /(x 2 )] = i
1 + 1 +
1
11
10'
The actual value is arctan2 « 1.107149, and our approximation is correct to two decimal places. H
8.2.1 How Good is the Composite Trapezoidal Rule?
We consider the composite trapezoidal rule for partitions of equal subintervals. Let Pi(x) be the
polynomial of degree < 1 that interpolates f(x) at Xj,Xj-|-i. Let
/i = /(*) dx, Tj = f i+1 Pi (x) dx = (xi+i - x.) ±^!±l) = £ {f{xi) +
That's right: the composite trapezoidal rule approximates the integral of fix) over [xj,Xj+i] by
the integral of pi(x) over the same interval.
Now recall our theorem on polynomial interpolation error. For x € [xj,Xj+i] , we have
f(x) - Pi{x) = — /( 2 )(£,.) (x - ^i) (a: - x i+1 ) ,
102
CHAPTER 8. INTEGRALS AND QUADRATURE
for some £ x G [x«, Xj+i] . Recall that ^ depends on x. To make things simpler, call it
Now integrate:
h-Ti= I f(x) - pi(x) dx = - / f"(€(x)) (x - Xi) (x - x i+1 ) dx.
J Xi ^ J Xi
We will now attack the integral on the right hand side. Recall the following theorem:
Theorem 8.7 (Mean Value Theorem for Integrals). Suppose / is continuous, g is Riemann
Integrable and does not change sign on [a, 0\. Then there is some ( G [a,0\ such that
f(x)g(x) dx = /(C) / g(x) dx.
We use this theorem on our integral. Note that (x — Xi) (x — Xj+i) is nonpositive on the interval
of question, [xj,Xj+i]. We assume continuity of f"(x), and wave our hands to get continuity of
f"(£(x)). Then we have
h-Ti = \f"{£) [ (x - Xi) {x - x i+1 ) dx,
Z J Xi
for some £j G By boring calculus and algebra, we find that
JXi
Xi+i /j3
(x - Xi) (x - x i+ i) dx = - — .
This gives
h 3
I'i ~ Ti = ~~^f
for some & G [xi,x i+1 ].
We now sum over all subintervals to find the total error of the composite trapezoidal rule
n— 1 , 3 ra— 1 /, \ 7 2
12 K ' 12
i=0 i=0
1 71— 1
-E/"(6
i=0
On the far right we have an average value, ^ Y17=o /"(£»)) which lies between the least and greatest
values of /" on the inteval [a, b], and thus by the IVT, there is some £ which takes this value. So
This gives us the theorem:
Theorem 8.8 (Error of the Composite Trapezoidal Rule). Let f"(x) be continuous on [a, b].
Let T be the value of the trapezoidal rule applied to f(x) on this interval with a partition of uniform
spacing, h, and let / = f(x) dx. Then there is some £ G [a, b] such that
Note that this theorem tells us not only the magnitude of the error, but the sign as well. Thus
if, for example, f(x) is concave up and thus /" is positive, then I — T will be negative, i.e., the
trapezoidal rule gives an overestimate of the integral /. See Figure 8.3.
8.2. TRAPEZOIDAL RULE
103
Figure 8.3: The trapezoidal rule is an overestimate for a function which is concave up, i.e., has
positive second derivative.
8.2.2 Using the Error Bound
Example Problem 8.9. How many intervals are required to approximate the integral
In 2 = I = / dx
Jo 1 + x
to within 1 x 10~ 10 ?
Solution: We have f(x) = thus f{x) = -j^- And f"{x) = Thus /"(£) is
continuous and bounded by 2 on [0, 1]. If we use n equal subintervals then Theorem 8.8 tells us the
error will be
i-o/i-o\ 2 no
12 \ n J 1 ^' 12n 2 '
To make this smaller than 1 x 10~ 10 , in absolute value, we need only take
i<lx 10- 10 ,
and so n > x 10 5 suffices. Because f"{x) is positive on this interval, the trapezoidal rule will
be an overestimate. H
Example Problem 8.10. How many intervals are required to approximate the integral
x 3 - 1 dx
o
to within 1 x 10 -6 ?
Solution: We have f(x) = x 3 — 1, thus f'{x) = 3x 2 , and f"(x) = 6x. Thus f"(£) is continuous
and bounded by 12 on [0, 2]. If we use n equal subintervals then by Theorem 8.8 the error will be
2-0 (2-Q\ 2 „ 2f"(0
12 V n ) 3n 2
104
CHAPTER 8. INTEGRALS AND QUADRATURE
To make this smaller than 1 x 10 6 , in absolute value, it suffices to take
24 «
s? 5 1 x 10
and so n > \/8 x 10 3 suffices. Because f"(x) is positive on this interval, the trapezoidal rule will
be an overestimate. H
8.3 Romberg Algorithm
Theorem 8.8 tells us, approximately, that the error of the composite trapezoidal rule approximation
is O (h 2 ). If we halve h, the error is quartered. Sometimes we want to do better than this. We'll
use the same trick that we did from Richardson extrapolation. In fact, the forms are exactly the
same.
Towards this end, suppose that /, a, b are given. For a given n, we are going to use the trapezoidal
rule on a partition of 2 n equal subintervals of [a, b}. That is h = ^r- Then define
lb - a 2 ™ _1
i=0
a
2 n
f(a) f(b) 2 ^\f .b-a
The intervals used to calculate 4>(n + 1) are half the size of those for (f)(n). As mentioned above,
this means the error is one quarter.
It turns out that if we had proved the error theorem differently, we would have proved the
relation:
rb
<K n ) = / f( x ) dx + a 2 h 2 l + a 4 /4 + a^h® + a%h % n + . . . ,
J a
where h n = ^pr- The constants are a function of f^\x) only. This should look just like something
from Chapter 7. What happens if we now calculate 4>(n + 1)? We have
rb
(f)(n + 1) = / f(x) dx + a 2 h 2 n+1 + a A h A n+1 + a 6 /i^ +1 + a 8 ^ +1 + . . . ,
J a
rb i 11 1
= J f(x) dx + -a 2 h 2 n + — a 4 h n + -^a G h % n + ^a 8 /i® + . . . .
This happens because h n+ \ = = = \- As with Richardon's method for approximating
derivatives, we now combine the right multiples of these:
/ 3 15 63
fix) dx + -4/4 + — a % h & n + — a 8 h 8 n + ...
/b 15 21
f{x) dx - -4/4 - r^a 6 h b n - ^a 8 hf
4o(n)-0(n+\) _ , , (i , ,,. ,
This approximation has a truncation error of O (/i^) .
Like in Richardson's method, we can use this to get better and better approximations to the
integral. We do this by constructing a triangular array of approximations, each entry depending
on two others. Towards this end, we let
R(n,0) = (f)(n),
8.3. ROMBERG ALGORITHM
105
then define, for m >
A m F?(n rn — 1\ — Ff(n — 1 m — 1\
(8.3)
, 4 m R(n, m - 1) - R(n - 1, m - 1)
R(n,m) =
4m _ 1
The familiar pyramid table then is:
12(0,0)
i?(l,0) i?(l,l)
R(2,0) R(2,l) R(2,2)
R(n,0) R(n,l) R(n,2) ■■■ R(n,n)
Even though this is exactly the same as Richardon's method, it has another name: this is called
the Romberg Algorithm.
Example Problem 8.11. Approximating the integral
dx
o 1 + ^
by Romberg's Algorithm; find R(l, 1).
Solution: The first column is calculated by the trapezoidal rule. Successive columns are found by
combining members of previous columns. So we first calculate R(0, 0) and R(l, 0). These are fairly
simple, the first is the trapezoidal rule on a single subinterval, the second is the trapezoidal rule on
two subintervals. Then
fl(0,0) = ^^[/(0) + /(2)] = jj,
R(l, 0) = 2 -^\ [/(0) + /(l) + /(l) + /(2)] =
Then, using Romberg's Algorithm we have
4fl(l,0)-i i (0,0) = t-j| = 32 =
v ' ; 4-1 3 30
H
At this point we are tempted to use Richardson's analysis. This would claim that R(n, n) is a
O (h^ n+1 ^ approximation to the integral. However, ho = b — a, and need not be smaller than 1.
This is a bit different from Richardson's method, where the original h is independently set before
starting the triangular array; for Romberg's algorithm, ho is determined by a and b.
We can easily deal with this problem by picking some k such that is small enough, say
smaller than 1. Then calculating the following array:
R(k,0)
R(k + 1,0) R{k + 1,1)
R{k + 2,0) J?(ife + 2,1) R(k + 2,2)
R(k + n,0) R(k + n,l) R(k + n,2) ••• R(k + n,n)
106
CHAPTER 8. INTEGRALS AND QUADRATURE
Quite often Romberg's Algorithm is used to compute columns of this array. Subtractive can-
celling or unbounded higher derivatives of f(x) can make successive approximations less accurate.
For this reason, entries in ever rightward columns are usually not calculated, rather lower entries
in a single column are calculated instead. That is, the user calculates the array:
R(k,0)
R(k + 1,0)
R(k + 2,0)
R(k + 1,1)
R(k + 2,1)
R(k + 2,2)
R(k + n,0) R(k + n,l) R(k + n,2)
R(k + n + 1,0) R(k + n + 1,1) R{k + n + 1,2)
R(k + n + 2,0) R(k + n + 2,1) R(k + n + 2,2)
R{k + n + 3,0) R(k + n + 3,l) R(k + n + 3,2)
R(k + n, n)
R{k + n + l,n)
R(k + n + 2,n)
R(k + n + 3,n)
Then R(k + n + l,n) makes a fine approximation to the integral as I — > oo. Usually n is small,
like 2 or 3.
8.3.1 Recursive Trapezoidal Rule
It turns out there is an efficient way of calculating R(n + 1, 0) given R(n, 0); first notice from the
above example that
R(Q,0) = t^Lh f{a) + mh
R(i,o) =
1 2
b-a l
2 2
f(a) + f(^) + f(^) + f(b)
It would be best to calculate R(l, 0) without recalculating f(a) and f{b). It turns out this is possible.
Let h n = and recall that
R{n, 0) = 4>{n) = h n
Thus
R(n + 1,0) = <t>(n + 1) = h n+1
/(«) + /(&)
2 n -l
+ Yl f( a + ih n)
i=l
2 n + l_ 1
1.
= -h,
2 + Z> f{a + ih n+1 )
/(«) + f(b)
2 n -l
+ J2f(a + (2i- l)h n+1 ) + / (a + (2i)h n ^j
/(«) + f(b)
2 n -l
2 n -l
+ f (a + */*n) + / (a + (2i - l)/i n+1 )
2 n -l
]-R{n, 0) + /i n+ i J] / (a + (2i - l)h n+l )
i=i
Then calculating + 1,0) requires only 2™ — 1 additional evaluations of f{x), instead of the
2 n+l + 1 usually required.
8.4. GAUSSIAN QUADRATURE
107
8.4 Gaussian Quadrature
The word quadrature refers to a method of approximating the integral of a function as the linear
combination of the function at certain points, i.e.,
/ f{x) dx « A f(x ) + A 1 f(x 1 ) + ... A n f{x n ),
J a
for some collection of nodes {xj}™ =0 , and weights {Ai}™ =Q . Normally one finds the nodes and weights
in a table somewhere; we expect a quadrature rule with more nodes to be more accurate in some
sense-the tradeoff is in the number of evaluations of /(■)■ We will examine how these rules are
created.
(8.4)
8.4.1 Determining Weights (Lagrange Polynomial Method)
Suppose that the nodes {xj}™ =0 are given. An easy way to find "good" weights {^4j}™ =0 for these
nodes is to rig them so the quadrature rule gives the integral of p(x), the polynomial of degree < n
which interpolates f[x) on these nodes. Recall
i=0
where £i(x) is the i th Lagrange polynomial. Thus our rigged approximation is the one that gives
rb rb n rb
j f(x)dx& / p(x) dx = y^J(xj) / li{x)dx.
i=0
If we let
r b
A.
= / £i{x)dx,
J a
then we have a quadrature rule.
If f(x) is a polynomial of degree < n then f(x) = p(x), and the quadrature rule is exact.
Example Problem 8.12. Construct a quadrature rule on the interval [0,4] using nodes 0, 1, 2.
Solution: The nodes are given, we determine the weights by constructing the Lagrange Polyno-
mials, and integrating them.
4>(S) = ^~^n"o' = ^ X ~ " 2 )'
'-..(-) - "rzr:: --ix)i,-2).
(x
-l)(x
-2)
(0
-1)(0
-2)
(x
- 0)(x
-2)
(1
-0)(1
-2)
(x
- 0)(x
-1)
(2
-0)(2
-1)
^(x) = ) n n ;; o ^ = -(x)(x - 1).
Then the weights are
A) = J to(x) dx = J i(x - l)(x - 2) dx = |,
f 4 f A 16
A 1 = / £ 1 (x)dx= / -(x)(x-2)dx = ,
Jo Jo 3
A 2 = £2(x)dx= / -(x)(x-l)dx
Jo Jo 2
20
108
CHAPTER 8. INTEGRALS AND QUADRATURE
Thus our quadrature rule is
f 8 16 20
^ /(x)dx«-/(0)- T /(l) + T /(2).
We expect this rule to be exact for a quadratic function f(x). To illustrate this, let f(x)
By calculus we have
x 2 + l.
Jo
x 2 + 1 dx= -x 3 + x
64 „ 76
y +4 = y
The approximation is
£, 2 + ld*^[0 + l]-^[l + l] + f [4 + 1] = f.
8.4.2 Determining Weights (Method of Undetermined Coefficients)
Using the Lagrange Polynomial Method to find the weights Ai is fine for a computer, but can
be tedious (and error-prone) if done by hand (say, on an exam). The method of undetermined
coefficients is a good alternative for finding the weights by hand, and for small n.
The idea behind the method is to find n + 1 equations involving the n + 1 weights. The equations
are derived by letting the quadrature rule be exact for f(x) = x J for j = 0, 1, . . . , n. That is, setting
rb n
/ x j dx = ^ AiixiY .
Ja i=0
For example, we reconsider Example Problem 8.12.
Example Problem 8.13. Construct a quadrature rule on the interval [0,4] using nodes 0, 1, 2.
Solution: The method of undetermined coefficients gives the equations:
J 1 dx = 4 = A + A 1 + A 2
Jo
4
x dx = 8 = A-y + 2A 2
4
x 2 dx = 64/3 = Ai +AA 2 .
We perform Naive Gaussian Elimination on the system:
1
1
(•
1
2
1
4
64/3 /
We get the same weights as in Example Problem 8.12: A 2 = ^-,Ai = — A) = §• ^
Notice the difference compared to the Lagrange Polynomial Method: undetermined coefficients
requires solution of a linear system, while the former method calculates the weights "directly."
Since we will not consider n to be very large, solving the linear system may not be too burdensome.
Moreover, the method of undetermined coefficients is useful in more general settings, as illus-
trated by the next example:
L
L
8.4. GAUSSIAN QUADRATURE
109
Example Problem 8.14. Determine a "quadrature" rule of the form
f 1 f(x)dx^Af(l) + Bf'(l) + Cf"(l)
Jo
that is exact for polynomials of highest possible degree. What is the highest degree polynomial for
which this rule is exact?
Solution: Since there are three unknown coefficients to be determined, we look for three equations.
We get these equations by plugging in successive polynomials. That is, we plug in f(x) = 1, x, x 2
and assuming the coefficients give equality:
f Idx = 1 = Al +B0+C0 = A
Jo
/ xdx = l/2 = Al + Bl + CO = A + B
Jo
[ x 2 dx = l/3 = Al + B2 + C2 = A + 2B + 2C
Jo
This is solved by A = 1, B = — 1/2, C = 1/6. This rule should be exact for polynomials of degree
no greater than 2, but it might be better. We should check:
\ x 3 dx = 1/4/ 1/2 = 1-3/2 + 1 = A1 + B3 + C6,
Jo
and thus the rule is not exact for cubic polynomials, or those of higher degree. H
8.4.3 Gaussian Nodes
It would seem this is the best we can do: using n + 1 nodes we can devise a quadrature rule that is
exact for polynomials of degree < n by choosing the weights correctly. It turns out that by choosing
the nodes in the right way, we can do far better. Gauss discovered that the right nodes to choose
are the n + 1 roots of the (nontrivial) polynomial, q(x), of degree n + 1 which has the property
f b
/ x k q(x)dx = (0<k<n).
J a
(If you view the integral as an inner product, you could say that q(x) is orthogonal to the polyno-
mials x k in the resultant inner product space, but that's just fancy talk.)
Suppose that we have such a q(x)-we will not prove existence or uniqueness. Let f(x) be a
polynomial of degree < 2n + 1. We write
f(x) =p(x)q(x) +r(x).
Both p(x),r(x) are of degree < n. Because of how we picked q(x) we have
J dx = 0.
J a
Thus
/ p(x)q(x)<
J a
f'O f'O f'O f'O
/ f(x) dx = p(x)q(x) dx + / r(x)dx = / r(x)dx.
J a J a J a J a
110
CHAPTER 8. INTEGRALS AND QUADRATURE
Now suppose that the Ai are chosen by Lagrange Polynomials so the quadrature rule on the
nodes Xi is exact for polynomials of degree < n. Then
n n n
^2 A if( X i) = ^2 A i [p{Xi)q(xi] + r(Xi)] = ^ A i r ( x i)-
i=0 i=0 i=0
The last equality holds because the x\ are the roots of q(x). Because of how the Ai are chosen we
then have
n n r-b f-b
S ^A i f{xi) = ^A i r{xi)=j r(x)dx = f(x)dx.
i=0 i=0 Ja Ja
Thus this rule is exact for f(x). We have (or rather, Gauss has) made quadrature twice as good.
Theorem 8.15 (Gaussian Quadrature Theorem). Let xi be the n + 1 roots of a (nontrivial)
polynomial, q(x), of degree n + 1 which has the property
rb
I x k q{x) dx = (0 < k < n) .
J a
Let Ai be the coefficients for these nodes chosen by integrating the Lagrange Polynomials. Then the
quadrature rule for this choice of nodes and coefficients is exact for polynomials of degree < 2n + 1.
8.4.4 Determining Gaussian Nodes
We can determine the Gaussian nodes in the same way we determine coefficients. The example is
illustrative
Example Problem 8.16. Find the two Gaussian nodes for a quadrature rule on the interval [0, 2].
Solution: We will find the function q(x) of degree 2, which is "orthogonal" to l,x under the inner
product of integration over [0,2]. Thus we let q(x) = cq + c\x + C2X 2 . The orthogonality condition
becomes
,2 ,2
/ lq(x) dx = xq(x) dx = that is,
Jo Jo
/ Co + c\x + C2X 2 dx = cqx + c\x 2 + C2X 3 dx =
Jo Jo
Evaluating these integrals gives the following system of linear equations:
8
2c + 2ci + -c 2 = 0,
g
2c + -ci + 4c 2 = 0.
This system is "under determined," that is, there are two equations, but three unknowns. Notice,
however, that if q(x) satisfies the orthogonality conditions, then so does q(x) = aq(x), for any real
number a. That is, we can pick the scaling of q(x) as we wish.
With great foresight, we "guess" that we want c 2 = 3. This reduces the equations to
2c + 2ci = -8,
2c + -ci = -12.
8.4. GAUSSIAN QUADRATURE
111
Simple Gaussian Elimination (cf. Chapter 3) yields the answer cq = 2, c\ = —6, C2 = 3.
Then our nodes are the roots of q(x) = 2 — Qx + 3x 2 . That is the roots
6 ± V36 - 24 ^3
6 ~ ~3~'
These nodes are a bit ugly. Rather than construct the Lagrange Polynomials, we will use the
method of undetermined coefficients. Remember, we want to construct Aq,A\ such that
L
2 f(x) dx « A f(l - ^) + A 1 f(l + ^)
is exact for polynomial /(x) of degree < 1. It suffices to make this approximation exact for the
"building blocks" of such polynomials, that is, for the functions 1 and x. That is, it suffices to find
A ,Ai such that
s:
L
1 dx = A + A 1
2 xdx = A (l-^-) + A 1 (l + ^-)
This gives the equations
2 = A + A x
2 = A)(l-^) + A 1 (l + ^)
This is solved by A§ = A\ = 1.
Thus our quadrature rule is
^/(x)dx«/(l-^) + /(l + ^
V3,
We expect this rule to be exact for cubic polynomials.
Example Problem 8.17. Verify the results of the previous example problem for f(x)
Solution: We have
o
2 2
f{x)dx = (1/4) x 4
4.
The quadrature rule gives
/(I- f ) + /(! + f)
3
+ 1
3
1 _j^ + 8 ?-?^U[l + 8^ + 8? + ^
3 9 27 3 9 27
= (2 - V3 - \/3/9) + (2 + + \/3/9) = 4
Thus the quadrature rule is exact for f(x).
112
CHAPTER 8. INTEGRALS AND QUADRATURE
8.4.5 Reinventing the Wheel
While it is good to know the theory, it doesn't make sense in practice to recompute these things all
the time. There are books full of quadrature rules; any good textbook will list a few. The simplest
ones are given in Table 8.1.
See also http : //mathworld . wolfram . com/Legendre-GaussQuadrature . html
n
Ai
x =
A = 2
1
Xq = -y/1/3
Aq = 1
Xl = a/173
A 1 = l
xq = -a/ 3 / 5
A = 5/9
2
xi =
A x = 8/9
x 2 = a/3/5
Ay = 5/9
Table 8.1: Gaussian Quadrature rules for the interval [—1,1]. Thus f^ 1 f(x)dx « Yl^o Af{xi),
with this relation holding exactly for all polynomials of degree no greater than 2n + 1.
Quadrature rules are normally given for the interval [—1,1]. On first consideration, it would
seem you need a different rule for each interval [a, b]. This is not the case, as the following example
problem illustrates:
Example Problem 8.18. Given a quadrature rule which is good on the interval [—1, 1], derive a
version of the rule to apply to the interval [a, b\.
Solution: Consider the substitution:
Then
b — a b + a b — a
x = 1-\ , so ax = at.
2 2 2
/*^ . . f 1 „ fb-a b + a\ b-a ,
Letting
2 2
if f(x) is a polynomial, g(t) is a polynomial of the same degree. Thus we can use the quadrature
rule for [—1, 1] on g(t) to evaluate the integral of f(x) over [a, b\. H
Example Problem 8.19. Derive the quadrature rules of Example Problem 8.16 by using the
technique of Example Problem 8.18 and the quadrature rules of Table 8.1.
Solution: We have a = 0, b = 2. Thus
b — a„(b — a b + a\ „, „.
rf*) = — /(—*+— )=/(* + !)■
To integrate f(x) over [0, 2], we integrate g(t) over [—1, 1]. The standard Gaussian Quadrature rule
approximates this as
1 g(t) dt « g(- a/T73) + g(V^) = /(I " + /(I + V 7 ^)-
-l
This is the same rule that was derived (with much more work) in Example Problem 8.16. H
8.4. GAUSSIAN QUADRATURE
113
Exercises
(8.1) Use the composite trapezoidal rule, by hand, to approximate
r 3
/ x 2 dx (= 9)
Jo
Use the partition {xi} 2 =0 = {0, 1,3} . Why is your approximation an overestimate?
(8.2) Use the composite trapezoidal rule, by hand, to approximate
1 i
X + 1
dx (= In 2 » 0.693)
Use the partition {xi} 3 =0 = {0, \, ^, l} . Why is your approximation an overestimate? (Check:
I think the answer is 0.7)
(8.3) Use the composite trapezoidal rule, by hand to approximate
L
1 A
dx.
l + x
2
Use n = 4 subintervals. How good is your answer?
(8.4) Use Theorem 8.8 to bound the error of the composite trapezoidal rule approximation of
J 2 x 3 dx with n = 10 intervals. You should find that the approximation is an overestimate.
(8.5) How many equal subintervals of [0, 1] are required to approximate Jq 1 cos x dx with error
smaller than 1 x 10~ 6 by the composite trapezoidal rule? (Use Theorem 8.8.)
(8.6) How many equal subintervals would be required to approximate
o dx.
o l + x 2
to within 0.0001 by the composite trapezoidal rule? (Hint: Use the fact that \f"(x)\ < 8 on
[0, 1] for f(x) = 4/(1 + x 2 ))
(8.7) How many equal subintervals of [2, 3] are required to approximate J* 2 3 e x dx with error smaller
than 1 x 10~ 3 by the composite trapezoidal rule?
(8.8) Simpson's Rule for quadrature is given as
/
J a
A x
f{x) dx « — [f(x ) + 4/0n) + 2f(x 2 ) + 4/(x 3 ) + • • • + 2/(x n _ 2 ) + 4/(x n _i) + f(x n )] ,
where Ax = (b — a)/n, and n is assumed to be even. Show that Simpson's Rule for n = 2 is
actually given by Romberg's Algorithm as R(l, 1). As such we expect Simpson's Rule to be
a O (h A ) approximation to the integral.
(8.9) Find a quadrature rule of the form
f 1 f(x) dx « Af(0) + Bf (1/2) + Cf(l)
Jo
that is exact for polynomials of highest possible degree. What is the highest degree polynomial
for which this rule is exact?
114
CHAPTER 8. INTEGRALS AND QUADRATURE
(8.10) Determine a "quadrature" rule of the form
/ f(x)dx^Af(0) + Bf'(-l) + Cf'(l)
J-i
that is exact for polynomials of highest possible degree. What is the highest degree polynomial
for which this rule is exact? (Since this rule uses derivatives of /, it does not exactly fit our
definition of a quadrature rule, but it may be applicable in some situations.)
(8.11) Determine a "quadrature" rule of the form
that is exact for polynomials of highest possible degree. What is the highest degree polynomial
for which this rule is exact?
(8.12) Consider the so-called order n Chebyshev Quadrature rule:
Find the weighting c n and nodes Xi for the case n = 2 and the case n = 3. For what order
polynomials are these rules exact?
(8.13) Find the Gaussian Quadrature rule with 2 nodes for the interval [1,5], i.e., find a rule
Before you solve the problem, consider the following questions: do you expect the nodes to
be the endpoints 1 and 5? do you expect the nodes to be arranged symmetrically around the
midpoint of the interval?
(8.14) Find the Gaussian Quadrature rule with 3 nodes for the interval [—1, 1], i.e., find a rule
To find the nodes xq, x\, X2 you will have to find the zeroes of a cubic equation, which could be
difficult. However, you may use the simplifying assumption that the nodes are symmetrically
placed in the interval [—1, 1].
(8.15) Write code to approximate the integral of a / on [a, b] by the composite trapezoidal rule on
n equal subintervals. Your m-file should have header line like:
function iappx = trapezoidal (f , a, b,n)
You may wish to use the code:
x = a .+ (b-a) .* (0:n) ./ n;
If f is defined to work on vectors element-wise, you can probably speed up your computation
by using
bigsum = 0.5 * ( f(x(l)) + f(x(n+l)) ) + sum( f(x(2:(n))) );
(8.16) Write code to implement the Gaussian Quadrature rule for n = 2 to integrate / on the
interval [a, b] . Your m-file should have header line like:
function iappx = gauss2 (f , a,b)
(8.17) Write code to implement composite Gaussian Quadrature based on code from the previous
problem. Something like the following probably works:
8.4. GAUSSIAN QUADRATURE
115
function iappx = gaussComp(f ,a,b,n)
% code to approximate integral of f over n equal subintervals of [a,b]
x = a .+ (b-a) .* (0:n) ./ n;
iappx = 0;
for i=l:n
iappx += gauss2(f ,x(i) ,x(i+l)) ;
end
Use your code to approximate the error function:
Compare your results with the octave/Matlab builtin function erf. (Try help erf in octave,
or see http : / /mathworld . wolfram . com/Erf . html)
The error function is used in probability. In particular, the probability that a normal random
variable is within z standard deviations from its mean is
Thus erf(l/\/2) ~ 0.683, and erf(2/\/2) ~ 0.955. These numbers should look familiar to you.
ed(z/y/2)
CHAPTER 8. INTEGRALS AND QUADRATURE
Chapter 9
Least Squares
9.1 Least Squares
Least squares is a general class of methods for fitting observed data to a theoretical model function.
In the general setting we are given a set of data
X
x
x 1
X n
y
Vo
yi
y n
and some class of functions, T. The goal then is to find the "best" / € T to fit the data to y = f(x).
Usually the class of functions T will be determined by some small number of parameters; the number
of parameters will be smaller (usually much smaller) than the number of data points. The theory
here will be concerned with defining "best," and examining methods for finding the "best" function
to fit the data.
9.1.1 The Definition of Ordinary Least Squares
First we consider the data to be two vectors of length n + 1. That is we let
x = [xq X\ ...
and y = [y y 1 ... y n ] T .
The error, or residual, of a given function with respect to this data is the vector r
That is
y - f(x).
r = [r n ... r n ] 1 , where n = y % - f{x{).
Our goal is to find / G T such that r is reasonably small. We measure the size of a vector
by the use of norms, which were explored in Subsection 1.4.1. The most useful norms are the £ p
norms. For a given p with < p < oo, the £ p norm of r is defined as
w ' = iS r, J
For the purposes of approximation, the easiest norm to use is the £ 2 norm:
Definition 9.1.1 ((Ordinary) Least Squares Best Approximant). The least-squares best
approximant to a set of data, x, y from a class of functions, T, is the function /* € T that
minimizes the £ 2 norm of the error. That is, if /* is the least squares best approximant, then
|y -/*(*)
mm \\y
/(*)
117
118
CHAPTER 9. LEAST SQUARES
We will generally assume uniqueness of the minimum. This method is sometimes called the ordinary
least squares method. It assumes there is no error in the measurement of the data x, and usually
admits a relatively straightforward solution.
9.1.2 Linear Least Squares
We illustrate this definition using the class of linear functions as an example, as this case is reason-
ably simple. We are assuming
T = = ax + b\ a,b £R) .
We can now think of our job as drawing the "best" line through the data.
By Definition 9.1.1, we want to find the function f(x) = ax + b that minimizes
\ 1/2
2 \
,i=0
This is a minimization problem over two variables, x and y. As you may recall from calculus class,
it suffices to minimize the sum rather than its square root. That is, it suffices to minimize
n
^ [axi + b- yi} 2 .
i=0
We minimize this function with calculus. We set
= — = ^2 2x k {ax k + b-y k )
a
k=0
0= lf = ^2 2 ( aX k + b -Vk)
k=0
These are called the normal equations. Believe it or not these are two equations in two unknowns.
They can be reduced to
^2 x l a + ^2 x kb = Yx k y k
^2x k a+ (n + l)b = ^yk
The solution is found by naive Gaussian Elimination, and is ugly. Let
'12
dn
= £4
= d 2 i
= Y, Xk
d22
= n+l
ei
= £ x k y k
e2
= £y fc
We want to solve
dna + d 12 b = e\
d 2 ia + d 2 2b = e 2
9.1. LEAST SQUARES
119
Gaussian Elimination produces
^22^11 — ^12^21
^22^11 — 6^12^21
The answer is not so enlightening as the means of finding the solution.
We should, for a moment, consider whether this is indeed the solution. Our calculations have
only shown an extrema at this choice of (a, 6); could it not be a maxima?
9.1.3 Least Squares from Basis Functions
In many, but not all cases, the class of functions, is the span of a small set of functions. This
case is simpler to explore and we consider it here. In this case J- can be viewed as a vector space
over the real numbers. That is, for f,g £ J 7 , and a, (3 G R, then af + (3g £ T , where the function
af is that function such that (af)(x) = af(x).
Now let {<7j(^)}™ =0 be a set of m + 1 linearly independent functions, i.e.,
cogo(x) + cigi(x) + . . . + c m g m (x) = Vx
c = Ci
Then we say that T is spanned by the functions {gj(x)}™ =Q if
f( x ) = ^2 c j9j(x) I Cj e R, j = 0, 1, . . . ,
m
> .
In this case the functions gj are basis functions for T . Note the basis functions need not be unique:
a given class of functions will usually have more than one choice of basis functions.
Example 9.1. The class T = {f(x) = ax + b \ a, b 6 R} is spanned by the two functions go(x) = 1,
and gi(x) = x. However, it is also spanned by the two functions go(x) = 2x + 1, and g± (x) = x — 1.
To find the least squares best approximant of T for a given set of data, we minimize the square
of the £ 2 norm of the error; that is we minimize the function
0(c o ,ci,...,c m ) = ^2
k=0
J2 c j9j{xk)
Vk
(9.1)
Again we set partials to zero and solve
£2
fc=0
Vk
9i{x k )
This can be rearranged to get
m
E
3=0
J2gj{xk)9i(x k )
cj = ^2yk9i{x k )
k=0
120
CHAPTER 9. LEAST SQUARES
If we now let
n n
di 'j = ^9j{x k )gi{x k ), ei = ^y fc ffi(x fe ),
k=0 k=0
Then we have reduced the problem to the linear system (again, called the normal equations):
doi
d()2 ■
dom
CO
eo
du
di2 ■
dim
Cl
ei
d-20
di\
d22 ■
d2m
C2
e-2
dmO
dmX
d m 2 ■
d mm
Cm
Cm
(9.2)
The choice of the basis functions can affect how easy it is to solve this system. We explore this
in Section 9.2. Note that we are talking about the basis {gjjjlo, and not exactly about the class
of functions T .
For example, consider what would happen if the system of normal equations were diagonal. In
this case, solving the system would be rather trivial.
Example Problem 9.2. Consider the case where m = 0, and go{x) = hix. Find the least squares
approximation of the data
X
0.50
0.75
1.0
1.50
2.0
2.25
2.75
3.0
y
-1.187098
-0.452472
-0.068077
0.713938
1.165234
1.436975
1.725919
1.841422
Solution: Essentially we are trying to find the c such that clnx best approximates the data. The
system of equation 9.2 reduces to the 1-D equation:
[S fc ln 2 x k ] c = Y, k y k lnx k
For our data, this reduces to:
4.0960c = 6.9844
so we find c = 1.7052. The data and the least squares interpolant are shown in Figure 9.1. H
Example 9.3. Consider the awfully chosen basis functions:
9o(x) = - l) x 2 - + 1
gi{x) = x 3 + (|-l) (x 2 + x) +1
where e is small, around machine precision.
Suppose the data are given at the nodes xq = 0, x\ = 1,X2 = — 1. We want to set up the normal
equations, so we compute some of the dij. First we have to evaluate the basis functions at the nodes
Xj. But this example was rigged to give:
g (x ) = 1, g {xi) = 0, g {x 2 ) = e
gi(x ) = 1, gi{x±) = e, gi{x 2 ) =
After much work we find we want to solve
1 + e 2 1
co
yo + m
1 1 + e 2 _
. c i _
_ yo + m _
9.2. ORTHONORMAL BASES
121
-2.5 1 1 1 1 1 1 1 1
0.5 1 1.5 2 2.5 3
Figure 9.1: The data of Example Problem 9.2 and the least squares interpolant are shown.
However, the computer would only find this if it had infinite precision. Since it does not, and since
e is rather small, the computer thinks e 2 = 0, and so tries to solve the system
" 1 1 "
CO
1 1
c l
. yo + m _
When yi ^ y2, this has no solution. Bummer.
This kind of thing is common in the method of least squares: the coefficients of the normal
equations include terms like
9i(xk)9j(x k ).
When the gi are small at the nodes x^, these coefficients can get really small, since we are squaring.
Now we draw a rough sketch of the basis functions. We find they do a pretty poor job of
discriminating around all the nodes.
9.2 Orthonormal Bases
In the previous section we saw that poor choice of basis vectors can lead to numerical problems.
Roughly speaking, if gi{xk) is small for some i's and /c's, then some dij can have a loss of precision
when two small quantities are multiplied together and rounded to zero.
Poor choice of basis vectors can also lead to numerical problems in solution of the normal
equations, which will be done by Gaussian Elimination.
Consider the case where T is the class of polynomials of degree no greater than m. For simplicity
we will assume that all Xi are in the interval [0, 1]. The most obvious choice of basis functions is
gj(x) = x 3 . This certainly gives a basis for J 7 , but actually a rather poor one. To see why, look
at the graph of the basis functions in Figure 9.3. The basis functions look too much alike on the
given interval.
A better basis for this problem is the set of Chebyshev Polynomials of the first kind, i.e.,
gj{x) = Tj(x), where
T (x) = 1, Ti(x) = x, T i+ i(x) = 2xTi(x) - Tj_i(x).
122
CHAPTER 9. LEAST SQUARES
Figure 9.2: The basis functions 50 0*0 = (f — l) x<1 ~ § x + lj an d 51 (z) = x 3 + (§ — l) (x 2 + x) + 1
are shown for e = 0.05. Note that around the three nodes 0, 1, —1, these two functions take nearly
identical values. This can lead to a system of normal equations with no solution.
Figure 9.3: The polynomials x l for i = 0, 1, ... ,6 are shown on [0, 1]. These polynomials make a
bad basis because they look so much alike, essentially.
These polynomials are illustrated in Figure 9.4.
9.2.1 Alternatives to Normal Equations
It turns out that the Normal Equations method isn't really so great. We consider other methods.
First, we define A as the n x m matrix defined by the entries:
aij = gj(xi), i = 0, 1,... ,n, j = 0,1,... ,m.
9.2. ORTHONORMAL BASES
123
Figure 9.4: The Chebyshev polynomials Tj(x) for j = 0,1, ... ,6 are shown on [0, 1]. These polyno-
mials make a better basis for least squares because they are orthogonal under some inner product.
Basically, they do not look like each other.
That is
" 3o Oo)
31(20)
32(2o) •
■ 9m (20)
30(21)
91 (21)
32(2l) •
• 9m (2l)
30(22)
31(22)
32(22) •
• 9rn (22)
30(23)
31(23)
32(23) •
• 9m (23)
30(24)
31(^4)
32(24) •
• 3m (24)
. 30 On)
91 (Xn)
32 (X n ) ■
9m\Xn)
We write it in this way since we are thinking of the case where n 3> m, so A is "tall."
After some inspection, we find that the Normal Equations can be written as:
A T Ac = A T y.
(9.3)
Now let the vector c be identified, in the natural way, with a function in J- '. That is c is
identified with f(x) = YlY=o c j9j( x )- Y° u should now convince yourself that
Ac = /(*).
And thus the residual, or error, of this function is r = y — Ac.
In our least squares theory we attempted to find that c that minimized
\\y-Ac\\l = {y- Ac) T (y- Ac)
We can see this as minimizing the Euclidian distance from y to A c. For this reason, we will have
that the residual is orthogonal to the column space of A, that is we want
A T r = A T (y- Ac) = 0.
124
CHAPTER 9. LEAST SQUARES
This is just the normal equations. We could rewrite this, however, in the following form: find c, r
such that
1 A "
r
y
A T
c
This is now a system of n + m variables and unknowns, which can be solved by specialized means.
This is known as the augmented form. We briefly mention that naive Gaussian Elimination is not
appropriate to solve the augmented form, as it turns out to be equivalent to using the normal
equations method.
9.2.2 Ordinary Least Squares in octave/Matlab
The ordinary least squares best approximant is so fundamental that it is hard-wired into oc-
tave/Matlab's system solve operator. The general equation we wish to solve is
Ac = y.
This system is overdetermined: there are more rows than columns in A, i.e., more equations than
unknown variables. However, in octave/Matlab, this is solved just like the case where A is square:
octave : 1> c = A \ y ;
Thus, unless specifically told otherwise, you should not implement the ordinary least squares solu-
tion method in octave/Matlab. Under the hood, octave/Matlab is not using the normal equations
method, and does not suffer the increased condition number of the normal equations method.
9.3 Orthogonal Least Squares
The method described in Section 9.1, sometimes referred to as "ordinary least squares," assumes
that measurement error is found entirely in the dependent variable, and that there is no error in
the independent variables.
For example, consider the case of two variables, x and y, which are thought to be related by an
equation of the form y = mx + b. A number of measurements are made, giving the two sequences
{xj}™ =1 and {yi}™ = i- Ordinary least squares assumes, that
yi = mxi + b + €i,
where ej, the error of the i th measurement, is a random variable. It is usually assumed that
E[ej] = 0, i.e., that the measurements are "unbiased." Under the further assumption that the
errors are independent and have the same variance, the ordinary least squares solution is a very
good one. 1
However, what if it were the case that
Hi = m(xi + &) + b + ei,
i.e., that there is actually error in the measurement of the x{! In this case, the orthogonal least
squares method is appropriate.
lr rhis means that the ordinary least squares solution gives an unbiased estimator of m and b, and, moreover, gives
the estimators with the lowest variance among all linear, a priori estimators. For more details on this, locate the
Gauss-Markov Theorem in a good statistics textbook.
9.3. ORTHOGONAL LEAST SQUARES
125
The difference between the two methods is illustrated in Figure 9.5. In Figure 9.5(a), the
ordinary least squares method is shown; it minimizes the sum of the squared lengths of vertical
lines from observations to a line, reflecting the assumption of no error in Xj. The orthogonal least
squares method will minimize the sum of the squared distances from observations to a line, as
shown in Figure 9.5(b).
i 1 1 1 1 < 1 1 i 1 1 1 1 < 1 1
1 2 3 4 5 6 1 2 3 4 5 6
(a) Ordinary Least Squares (b) Orthogonal Least Squares
Figure 9.5: The ordinary least squares method, as shown in (a), minimizes the sum of the squared
lengths of the vertical lines from the observed points to the regression line. The orthogonal least
squares, shown in (b), minimizes the sum of the squared distances from the points to the line. The
offsets from points to the regression line in (b) may not look orthogonal due to improper aspect
ratio of the figure.
We construct the solution geometrically. Let {x^}^ be the set of observations, expressed as
vectors in ]R n . The problem is to find the vector n and number d such that the hyperplane n T x = d
is the plane that minimizes the sum of the squared distances from the points to the plane. We can
solve this problem using vector calculus.
The function
1 m n 2
J-Y\\n T x^-d .
is the one to be minimized with respect to n and d. It gives the sum of the squared distances from
the points to the plane described by n and d. To simplify its computation, we will minimize it
1 1 1 2
subject to the constraint that ||n|| 2 = 1. Thus our problem is to find n and d that solve
min n T £C^ — d
We can express this as min f(n,d) subject to g(n,d) = 1. The solution of this problem uses the
Lagrange multipliers technique.
Let £(n,d,\) = f(n,d) — Xg(n,d). The Lagrange multipliers theory tells us that a necessary
126
CHAPTER 9. LEAST SQUARES
condition for n, d to be the solution is that there exists A such that the following hold:
(§§(n,d,A)=0
V„£ (n, d, A) =
K g(n,d) = 1
Let X be the m x n matrix whose i th row is the vector »W . We can rewrite the Lagrange
function as
£ (n, cZ, A) = (Xn — (il m ) (Xn — dl m ) — An T n
= n T X T Xn - 2dl™ T Xn + d 2 l m T l m - Xn T n
We solve the necessary condition on
= — (n, d, A) = 2dl m T l m - 21 m T Xn
Xn/1
rn Irr?,
This essentially tells us that we will zero the first moment of the data around the line. Note that
this determination of d requires knowledge of n. However, since this is a necessary condition, we
can plug it into the gradient equation:
= V n £ (n, d, A) = 2X T Xn - 2d (l m T x) T - 2An = 2 X T Xn - hlL^lH (
X
An
The middle term is a scalar times a vector, so the multiplication can be commuted, this gives
X T Xn
X T
1 1 1
-n — Xn
vTi -i Tv
X T X - m } m - Al
1 1 1
TJ,
Thus n is an eigenvector of the matrix
M=X T X- XTl ™ r 1 ™ TX
1 T 1
J-m J -m
The final condition of the three Lagrange conditions is that g{n,d) = n n = 1. Thus if n is a
minimizer, then it is a unit eigenvector of M.
It is not clear which eigenvector is the right one, so we might have to check all the eigenvectors.
However, further work tells us which eigenvector it is. First we rewrite the function to be minimized,
when the optimal d is used: f(n) = f(n, l TTl T Xn/l TTl T l rn ,). We have
2
f(n) = n T X T Xn - 2
Xn.
1 T l
J-m J -m
Xn +
lm T Xn
1 T l
J-m J -m
1 T l
Tv Tv « T X T l m l m T Xn n T X T l m l m T Xnl m T l,
n X Xn — 2 : — — h -
n T X T Xn - 2
n T X T Xn
1 T l
n T X T l rT7 l rT7 T Xn n
+
1 T l 1 T l
T X T l m l m T Xn
1 T l
n T X T l m l m T Xn
n
T
x'x
X T l m l m T X
I Ti
n
n T Mn
9.3. ORTHOGONAL LEAST SQUARES
127
This sequence of equations only required that the d be the optimal d for the given n, and used the
fact that a scalar is it's own transpose, thus l m T Xn = ra T X T l m .
Now if n is a unit eigenvector of M, with corresponding eigenvalue A, then n T Mn = A. Note
that because f(n,d) is defined as w T w, for some vector w, then the matrix M must be positive
definite, i.e., A must be positive. However, we want to minimize A. Thus we take n to be the unit
eigenvector associated with the smallest eigenvalue.
Note that, by roundoff, you may compute that M has some negative eigenvalues, though they
are small in magnitude. Thus one should select, as n, the unit eigenvector associated with the
smallest eigenvalue in absolute value.
Example Problem 9.4. Find the equation of the line which best approximates the 2D data
{(xi,yi)}™ =1 , by the orthogonal least squares method.
Solution: In this case the matrix X is
Xl
Vi
X2
2/2
^3
2/3
Thus we have that
m = x t x _x^o^x
1
m
— m
x xy
xy y 2
E x i E Xiy,
E x iVt E vl
E x i E x iV-
E x iVi E yf
E x\ — rax 1 E x iy-i — mxy
E Xiyi - mxy E yf ~ m y 2
2a (3
13 2 7
(E Xi) E Xi E yi
E Xi E yi (E yiY
(9.4)
where we use x, y, to mean, respectively, the mean of the x- and y-values. The characteristic
polynomial of the matrix, whose roots are the eigenvalues of the matrix is
p(A) = det
2a- A
f3 2 7 -A
= 4a 7 -/3 2 -2(a + 7 )A + A 2 .
The roots of this polynomial are given by the quadratic equation
2(q + 7) ± ^4(a + 7 ) 2 -4(4a 7 -/3 2 )
A±
(a + 7)± V(«-7) +/3 2
We want the smaller eigenvalue, A_ = (a + 7) — ^/ (a — r y) 2 + j3 2 . The associated eigenvector is in
the null space of the matrix
=
2a - A_
P 2 7 - A_
v =
2a - A_ /3
P 2 7 -A.
-A:
1
/?- k(2a- A_)
-A;/? + 27- A_
128
CHAPTER 9. LEAST SQUARES
This is solved by
2 7 - A_ _ (7 " ") + \Jh~af + P 2
k =
(9.5)
YsiVi -xj)-m (y 2 - x 2 ) + y {Y,{yf ~ tf) -m(y 2 - x 2 )) 2 + 4 (]T Xi yi - mxy) 2
2E Xiyi ~ mxy
Thus the best plane through the data is of the form
— kx + ly = d or y = kx + d
Note that the eigenvector, v, that we used is not a unit vector. However it is parallel to the unit
vector we want, and can still use the regular formula to compute the optimal d.
d — l m T Xn/l T7T T l T7T — — lm. T X
m
-k
1
[ x y
-k
1
y — kx
Thus the optimal line through the data is
y = kx + d= k(x — x)+y or y — y = k(x — x) ,
with k defined in equation 9.5.
9.3.1 Computing the Orthogonal Least Squares Approximant
Just as the Normal Equations equation 9.3 define the solution to the ordinary least squares ap-
proximant, but are not normally used to solve the problem, so too is the above described means of
computing the orthogonal least squares not a good idea.
First we define
W = X
1 7n 1 rrv, X
I Ti •
Now note that
W W - i X - lml ™ TX ) (x - lml r TX ) X X - 2
1 1 1
1 1 1
X T lmlm T X (lmlm^X) l^l^^X
lm T lrTi (l m T lm) 2
V T V r) X T l T 7alTTi T X X T l Tn ,(l rn T l rn ,)l TT1 T X j X T l m l m T X X T l Tn ,l T17 , T X
— A A — Z =p 1 — ^ — A A — z =p 1 =p
1 T l
(1 T i V
1 T l
J-tti - L m.
1 T l
Y T 1 1 Tv
= X T X- ™ m = M.
1 T 1
We now need some linear algebra magic:
Definition 9.3.1. A square matrix U is called unitary if its inverse is its transpose, i.e.,
U T U = I = UU T .
From the first equation, we see that the columns of U form a collection of orthonormal vectors.
The second equation tells that the rows of U are also such a collection.
9.3. ORTHOGONAL LEAST SQUARES
129
Definition 9.3.2. Every m x n matrix B can be decomposed as
B = UZV T ,
where U and V are unitary matrices, U is m x m, V is n x n, and Y. is m x n, and has nonzero
elements only on its diagonal. The values of the diagonal of T. are the singular values of B. The
column vectors of V span the row space of B, and the column vectors of U contain the column space
of B.
The singular value decomposition generalizes the notion of eigenvalues and eigenvalues to the
case of nonsquare matrices. Let be the i th column vector of U, the i th column of V, and
let <7j be the i th element of the diagonal of Z. Then if x = J2i otiv( l \ we have
Bx = UZV
• CXr.
G\OL\
O20L2
OiOLiU
(0
Thus t>W and act as an "eigenvector pair" with "eigenvalue" cij.
The singular value decomposition is computed by the octave/Matlab command [U,S,V] =
svd(B). The decomposition is computed such that the diagonal elements of S, the singular values,
are positive, and decreasing with increasing index.
Now we can use the singular value decomposition on W:
W = UZV
T
Thus
M = W T W = VZ T U T UZV T = VZ T IZV T = VS 2 V T ,
where S 2 is the n x n diagonal matrix whose i th diagonal is of. Thus we see that the solution to
the orthogonal least squares problem, the eigenvector associated with the smallest eigenvalue of M,
is the column vector of V associated with the smallest, in absolute value, singular value of W. In
octave/Matlab, this is V(: ,n).
This may not appear to be a computational savings, but if M is computed up front, and its
eigenvectors are computed, there is loss of precision in its smallest eigenvalues which might lead us
to choose the wrong normal direction (especially if there is more than one!). On the other hand,
M is a n x n matrix, whereas W is m x n. If storage is at a premium, and the method is being
reapplied with new observations, it may be preferrable to keep M, rather than keeping W. That is,
it may be necessary to trade conditioning for space.
9.3.2 Principal Component Analysis
The above analysis leads us to the topic of Principal Component Analysis. Suppose the m x n
matrix X represents m noisy observations of some process, where each observation is a vector in
W 1 . Suppose the observations nearly lie in some /c-dimensional subspace of R n , for k < n. How can
you find an approximate value of k, and how do you find the /c-dimensional subspace?
130
CHAPTER 9. LEAST SQUARES
As in the previous subsection, let
W = X - hn^L^ = X - l m X. where X = l m T X/l m T l m .
Now note that X is the mean of the m observations, as a row vector in M n . Under the assumption
that the noise is approximately the same for each observation, we should think that X is likely
to be in or very close to the fc-dimensional subspace. Then each row of W should be like a vector
which is contained in the subspace. If we take some vector v and multiply Wd, we expect the
resultant vector to have small norm if v is not in the /c-dimensional subspace, and to have a large
norm if it is in the subspace.
You should now convince yourself that this is related to the singular value decomposition of W.
Decomposing v = Y2i a i v ^\ we have \Nv = dionu^, and thus product vector has small norm if
the a.% associated with large cij are small, and large norm otherwise. That is the principal directions
of the "best" fc-dimensional subspace associated with X are the k columns of V associated with the
largest singular values of W.
0.3 1 1 1 1 1 1 1 1 1
5 10 15 20 25 30 35 40
Figure 9.6: Noisy 5-dimensional data lurking in M are treated to Principal Component Analysis.
This figure shows the normalized cumulative sum of squared singular values of W. There is a sharp
"elbow" at k = 5, which indicates the data is actually 5-dimensional.
If k is unknown a priori, a good value can be obtained by eyeing a graph of the cumulative
sum of squared singular values of W, and selecting the k that makes this appropriately large. This
is illustrated in Figure 9.6, with data generated by the following octave/Matlab code:
octave :1> X = rand(150,5) ;
octave:2> sawX = (X * randn(5,40)) + kron(ones(150, 1) ,randn(l,40)) ;
octave :3> sawX += 0.1 * randn(150,40) ;
octave:4> wun = ones(150,l);
9.3. ORTHOGONAL LEAST SQUARES
131
octave :5> W = sawX - (wun * wun' * sawX) / (wun' * wun) ;
octave:6> [U,S,V] = svd(W) ;
octave :7> eigs = diag(S'*S);
octave:8> cseig = cumsum(eigs) ./ sum(eigs) ;
octave :9> plot ( (1 : 40) ' , cseig, 'xr-@ ; cum. sum sqrd. sing, vals.;');
132
CHAPTER 9. LEAST SQUARES
Exercises
(9.1) How do you know that the choice of constants Cj in our least squares analysis actually find a
minimum of equation 9.1, and not, say, a maximum?
(9.2) Our "linear" least squares might be better called the "affine" least squares. In this exercise
you will find the best linear function which approximates a set of data. That is, find the
function f(x) = cx which is the ordinary least squares best approximant to the given data
X
x
Xi
y
yo
Vi
Vn
(9.3) Find the constant that best approximates, in the ordinary least squares sense, the given data
x, y. (Hint: you can use equation 9.2 or equation 9.3 using a single basis function go(x) = 1.)
Do the X{ values affect your answer?
(9.4) Find the function f(x) = c which best approximates, in the ordinary least squares sense, the
data
X
l
-2
5
y
l
-2
4
(9.5) Find the function ax + b that best approximates the data
X
-l
2
y
l
-l
Use the ordinary least squares method.
(9.6) Find the function ax + b that best approximates the data of the previous problem using the
orthogonal least squares method.
(9.7) Find the function ax 2 + b that best approximates the data
X
1
2
y
0.3
0.1
0.5
(9.8) Prove that the least squares solution is the solution of the normal equations, equation 9.3,
and thus takes the form
c= ^A T A^ 1 A T y.
(9.9) For ordinary least squares regression to the line y = mx + b, express the approximate m
in terms of the a, j3, and 7 parameters from equation 9.4. Compare this to that found in
equation 9.5.
(9.10) Any symmetric positive definite matrix, G can be used to define a norm in the following
way:
\\v\\ G =df V v T Gv
The weighted least squares method for G and data A and y, finds the c that solves
min || Ac — y||^
Prove that the normal equations form of the solution of this problem is
A T GAc = A T Gy.
9.3. ORTHOGONAL LEAST SQUARES
133
(9.11) (Continuation) Let the symmetric positive definite matrix G have Cholesky Factorization
G = LL T , where L is a lower triangular matrix. Show that the solution to the weighted least
squares method is the c that is the ordinary (unweighted) least squares best approximate
solution to the problem
L T Ac = L T i
(9.12) The r 2 statistic is often used to describe the "goodness-of-fit" of the ordinary least squares
approximation to data. It is defined as
r
2 l|Ac|| 2
II l|2 '
Hi/112
where c is the approximant A T A 1 A T y.
(a) Prove that we also have
r
1 _ 1 _ l|Ac — y|||
II l|2
(b) Prove that the r 2 parameter is "scale- invariant," that is, if we change the units of our
data, the parameter is unchanged (although the value of c may change.)
(c) Devise a similar statistic for the orthogonal least squares approximation which is scale-
invariant and rotation-invariant.
(9.13) As proved in Example Problem 9.4, the orthogonal least squares best approximate line for
data {{xi,yi)}™ =1 goes through the point (x,y). Does the same thing hold for the ordinary
least squares best approximate line?
(9.14) Find the constant c such that f(x) = In (cx) best approximates, in the least squares sense,
the given data
X
xo
X\
Xn
y
yo
yi
Vn
(Hint: You cannot use basis functions and equation 9.2 to solve this. You must use the
Definition 9.1.1.) The geometric mean of the numbers a\, 0,2, . . . ,a n is defined as (Y\ a-i) l ^ n ■
How does your answer relate to the geometric mean?
(9.15) Write code to find the function ae x + be~ x that best interpolates, in the least squares sense,
a set of data {(xi, yi)Y~^. Your m-file should have header line like:
function [a,b] = lsqrexp(xys)
where xys is the (n + 1) x 2 matrix whose rows are the n + 1 tuples (xi,yi). Use the normal
equations method. This should reduce the problem to solving a 2 x 2 linear system; let
octave/Matlab solve that linear system for you. (Hint: To find x such that Mx = b, use x =
M \ b.)
(a) What do you get when you try the following?
octave :1> xs = [-1 -0.5 0.5 1]';
octave : 2> ys = [1.194 0.430 0.103 0.322 1.034]';
octave :3> xys = [xs ys] ;
octave :4> [a,b] = lsqrexp(xys)
(b) Try the following:
octave :5> xys = xys = [-0.00001 0.3;0 0.6; 0.00001 0.7];
octave :6> [a,b] = lsqrexp(xys)
Does something unexpected happen? Why?
(9.16) Write code to find the plane n T x^ = d that best interpolates, in the least squares sense, a
set of data {x^Y i= ™- Your m-file should have header line like:
CHAPTER 9. LEAST SQUARES
function [n,d] = orthlsq(X)
where X is the m x n matrix whose rows are the m tuples (xi, X2, ■ ■ • , x n ). Use the oc-
tave/Matlab eig function, which returns the eigenvalues and vectors of a matrix.
(a) Create artificially planar data by the following:
octave:l> m = 100;n = 5;epsi = 0.1;offs = 2;
octave :2> Xsub = rand(m, (n-1) ) ;
octave :3> T = rand (n-1 ,n) ;
octave :4> X = Xsub * T;
octave :5> X += offs .* ones (size (X) ) ;
octave :6> X += epsi .* randn(size(X) ) ;
Now X contains data which are approximately on a (n — l)-dimensional hyperplane in
W l . The data, when centered, lay in the rowspace of T. Now test your code.
octave:7> [n,d] = orthlsq(X) ;
You should have that the vector n is orthogonal to the rows of T. Check this:
octave :8> disp(norm(T * n));
What do you get? Try again with epsi very small (1 x 10 ~ 5 ) or zero.
(b) Compare the orthogonal least squares technique to the ordinary least squares in a two
dimensional setting. Consider x, y data nearly on the line y = mx + b:
octave :1> obs = 100;xeps = 1.0;yeps = 1.0 ;m = 0.375;b = 1;
octave :2> Xtrue = randn(obs , 1) ; Ytrue = m .* Xtrue .+ b;
octave :3> XSaw = Xtrue + xeps * randn(size (Xtrue) ) ;
octave :4> YSaw = Ytrue + yeps * randn (size (Ytrue) ) ;
octave :5> lsmb = [XSaw, ones (obs, 1)] \ YSaw;
octave :6> Ism = lsmb(l); lsb = lsmb (2) ;
octave :7> [n,d] = orthlsq( [XSaw, YSaw]);
octave :8> orthm = -n(l) / n(2);orthb = d / n(2) ;
Compare the estimated coefficients from both solutions. Plot the two regression lines,
octave :9> hold off;
octave: 10> plot (XSaw, YSaw, "*; observed data;");
octave :11> hold on;
octave: 12> plot (XSaw, m*XSaw .+ b,"r-;true line;");
octave :13> plot (XSaw, lsm*XSaw .+ lsb, "g-; ordinary least squares;");
octave :14> plot (XSaw, orthm*XSaw .+ orthb, "b-; orthogonal least squares;");
octave :15> hold off;
Try this again with different values of xeps and yeps, which control the sizes of the
errors in the x and y dimensions. Try with xeps = 0.1, yeps = 1.0, and with with
xeps = 1.0, yeps = 0.1. Under which situations does orthogonal least squares outperform
ordinary least squares? Do they give equally erroneous results in some situations? Can
you think of a way of "reweighting" channel data to make the orthogonal least squares
method more robust in cases of unequal error variance?
Chapter 10
Ordinary Differential Equations
Ordinary Differential Equations, or ODEs, are used in approximating some physical systems. Some
classical uses include simulation of the growth of populations and trajectory of a particle. They are
usually easier to solve or approximate than the more general Partial Differential Equations, PDEs,
which can contain partial derivatives.
Much of the background material of this chapter should be familiar to the student of calculus.
We will focus more on the approximation of solutions, than on analytical solutions. For more
background on ODEs, see any of the standard texts, e.g., [8].
10.1 Elementary Methods
A one-dimensional ODE is posed as follows: find some function x(t) such that
In calculus classes, you probably studied the case where f(t,x(t)) is independent of its second
variable. For example the ODE given by
A more general case is when f(t,x) is separable, that is f(t,x) = g(t)h{x) for some functions
g, h. You may recall that the solution to such a separable ODE usually involved the following
equation:
Finding an analytic solution can be considerably more difficult in the general case, thus we turn
to approximation schemes.
(10.1)
has solution x(t) = — ^t 2 + K, where K is chosen such that x(a) = c.
135
136
CHAPTER 10. ORDINARY DIFFERENTIAL EQUATIONS
10.1.1 Integration and 'Stepping'
We attempt to solve the ODE by integrating both sides. That is
*^± = f(t,x(t)), yields
rt+h r t+h
/t+n rt+n
dx = f(r,x(r))dr, thus
/t+h
f{r,x{r))dr. (10.2)
If we can approximate the integral then we have a way of 'stepping' from t to t + h, i.e., if we have
a good approximate of x(t) we can approximate x(t + h). Note that this means that all the work
we have put into approximating integrals can be used to approximate the solution of ODEs.
Using stepping schemes we can approximate x(tfi na i) given x(ti n i t i a i) by taking a number of
steps. For simplicity we usually use the same step size, h, for each step, though nothing precludes
us from varying the step size.
If we apply the left-hand rectangle rule for approximating integrals to equation 10.2, we get
x(t + h)^x(t) + hf(t,x(t)). (10.3)
This is called Euler's Method. Note this is essentially the same as the "forward difference" approx-
imation we made to the derivative, equation 7.1.
Trapezoid rule gives
x(t + /i) » x{t) + ^ [f(t, x{t)) + f(t + h, x(t + h))\ .
But we cannot evaluate this exactly, since x(t + h) appears on both sides of the equation. And it is
embedded on the right hand side. Bummer. But if we could approximate it, say by using Euler's
Method, maybe the formula would work. This is the idea behind the Runge-Kutta Method (see
Section 10.2).
10.1.2 Taylor's Series Methods
We see that our more accurate integral approximations will be useless since they require information
we do not know, i.e., evaluations of f(t,x) for yet unknown x values. Thus we fall back on
Taylor's Theorem (Theorem 1.1). We can also view this as using integral approximations where all
information comes from the left-hand endpoint.
By Taylor's theorem, if x has at least m + 1 continuous derivatives on the interval in question,
we can write
x (t + h)= x{t) + hx'{t) + \h 2 x"{t) + ... + —h m x^ m \t) + - — l —-h m+1 x^ m+1 \T),
where r is between t and t + h. Since r is essentially unknown, our best calculable approximation
to this expression is
x{t + h)r* x(t) + hx'{t) + \h 2 x"{t) + ... + —h m x^ (t).
This approximate solution to the ODE is called a Taylor's series method of order m. We also say
that this approximation is a truncation of the Taylor's series. The difference between the actual
10.1. ELEMENTARY METHODS
137
value of x(t + h), and this approximation is called the truncation error. The truncation error exists
independently from any error in computing the approximation, called roundoff error. We discuss
this more in Subsection 10.1.5.
10.1.3 Euler's Method
When m = 1 we recover Euler's Method:
x{t + h) = x(t) + hx'(t) = x{t) + hf(t, x{t)).
Example 10.1. Consider the ODE
dt ~ x '
z(0) = 1.
The actual solution is x{t) = e*. Euler's Method will underestimate x(t) because the curvature of
the actual x(t) is positive, and thus the function is always above its linearization. In Figure 10.1,
we see that eventually the Euler's Method approximation is very poor.
fculer approxir^atioij
1 1 1 1 1 1 1
0.5 1 1.5 2 2.5 3
Figure 10.1: Euler's Method applied to approximate x' = x, x(0) = 1. Each step proceeds on
the linearization of the curve, and is thus an underestimate, as the actual solution, x(t) = e has
positive curvature. Here step size h = 0.3 is used.
10.1.4 Higher Order Methods
Higher order methods are derived by taking more terms of the Taylor's series expansion. Note
however, that they will require information about the higher derivatives of x(t), and thus may not
be appropriate for the case where f(t,x) is given as a "black box" function.
138
CHAPTER 10. ORDINARY DIFFERENTIAL EQUATIONS
Example Problem 10.2. Derive the Taylor's series method for m = 3 for the ODE
dx( fil — T I p x
At '
x(0) = 1.
Solution: By simple calculus:
x'(t) = x + e x
x "(t) = x' + e x x'
x"'(t) = x"(l + e x ) + x'e x x'
We can then write our step like a program:
x'(t) <- x(t) + e x(t)
z"(t) «-x'(t) (l + e x ^
x »>(t) <_ x "( t ) (l + e -W) + V«
x(/ + /i) « x(t) + hx'(t) + J/iV(t) + -/i 3 x'"(t)
— ■ D
10.1.5 Errors
There are two types of errors: truncation and roundoff. Truncation error is the error of truncating
the Taylor's series after the m th term. You can see that the truncation error is O (/i m+1 ). Roundoff
error occurs because of the limited precision of computer arithmetic.
While it is possible these two kinds of error could cancel each other, given our pessimistic
worldview, we usually assume they are additive. We also expect some tradeoff between these two
errors as h varies. As h is made smaller, the truncation error presumably decreases, while the
roundoff error increases. This results in an optimal /i-value for a given method. See Figure 10.2.
Unfortunately, the optimal /i-value will usually depend on the given ODE and the approximation
method. For an ODE with unknown solution, it is generally impossible to predict the optimal h. We
have already observed this kind of behaviour, when analyzing approximate derivatives-c/. Example
Problem 7.5.
Both of these kinds of error can accumulate: Since we step from x(t) to x(t + h), an error in
the value of x(t) can cause a greater error in the value of x(t + h). It turns out that for some ODEs
and some methods, this does not happen. This leads to the topic of stability.
10.1.6 Stability
Suppose you had an exact method of solving an ODE, but had the wrong data. That is, you were
trying to solve
x(a) = c,
but instead fed the wrong data into the computer, which then solved exactly the ODE
^ =/&*(*))>
x(a) = c + e.
10.1. ELEMENTARY METHODS
139
truncation _ _ _ _
roundoff
1e-07 1e-06 1e-05 0.0001 0.001 0.01
Figure 10.2: The truncation, roundoff, and total error for a fictional ODE and approximation
method are shown. The total error is apparently minimized for an h value of approximately 0.0002.
Note that the truncation error appears to be O (h), thus the approximation could be Euler's Method.
In reality, we expect the roundoff error to be more "bumpy," as seen in Figure 7.1.
This solution can diverge from the correct one because of the different starting conditions. We
illustrate this with the usual example.
Example 10.3. Consider the ODE
dx(t)
— ; — = x -
dt
This has solution x(t) = x(0)e . The curves for different starting conditions diverge, as in Fig-
ure 10.3(a).
However, if we instead consider the ODE
dx(t)
which has solution x(t) = x(0)e , we see that differences in the initial conditions become immate-
rial, as in Figure 10.3(b).
Thus the latter ODE exhibits stability: roundoff and truncation errors accrued at a given step
will become irrelevant as more steps are taken. The former ODE exhibits the opposite behavior-
accrued errors will be amplified.
140
CHAPTER 10. ORDINARY DIFFERENTIAL EQUATIONS
Figure 10.3: Exponential functions with different starting conditions: in (a), the functions (l + e)e
are shown, while in (b), (l+e)e~* are shown. The latter exhibits stability, while the former exhibits
instability.
10.1. ELEMENTARY METHODS
141
It turns out there is a simple test for stability. If f x > 6 for some positive 8, then the ODE
is instable; if f x < —5 for a positive 5, the equation is stable. Some ODEs fall through this test,
however.
We can prove this without too much pain. Define x(t, s) to be the solution to
^ = f(t,x(t)),
x(a) = s,
for t > a.
Then instability means that
d
lim
t^oo
oo.
Think about this in terms of the curves for our simple example x' = x.
If we take the derivative of the ODE we get
d
— x(t,s) = f(t,x(t,s)),
^| x(M) = ^ /(t ' x(M)) ' Z ' e "
^"^U(t, S) = f t (t, S(t))^ + f x (t, X(t, S)) 9X{
os at as as
By the independence of t, s the first part of the RHS is zero, so
||x(M) = / ,(t,x(M))^M.
We use continuity of x to switch the orders of differentiation:
9 9, \ , , , u 9 , .
— — x{t,s) = f x (t,x{t,s))—x{t,s).
Defining u(t) = -^x(t,s), and q(t) = f x (t,x(t,s)), we have
^u(t)=q(t)u(t).
The solution is u(t) = ce^\ where
Q(t)= f q{r)dr.
J a
If q(r) = f x (r,x(r,s)) > 5 > 0, then lim Q(t) = oo, and so limu(i) = oo. But note that u(t) =
■§^x(t, s), and thus we have instability.
Similarly if q(r) < — 5 < 0, then limQ(t) = — oo, and so limu(t) = 0, giving stability.
10.1.7 Backwards Euler's Method
Euler's Method is the most basic numerical technique for solving ODEs. However, it may suffer
from instability, which may make the method inappropriate or impractical. However, it can be
recast into a method which may have superior stability at the cost of limited applicability.
142
CHAPTER 10. ORDINARY DIFFERENTIAL EQUATIONS
This method, the so-called Backwards Euler's Method, is a stepping method: given a reasonable
approximation to x(t), it calculates an approximation to x(t + h). As usual, Taylor's Theorem is
the starting point. Letting tf = t + h, Taylor's Theorem states that
x(t f -h) = x(t f ) + (-h)x'(t f ) + ^-x"(t f ) + ...,
for an analytic function. Assuming that x has two continuous derivatives on the interval in question,
the second order term is truncated to give
x(t) =x{t + h-h)= x(t f - /i) ~ x(t f ) + (-h)x'(t } )
and so, x(t + h) « x(t) + hx'(t + h)
The complication with applying this method is that x'(t+h) may not be computable if x(t+h) is
not known, thus cannot be used to calculate an approximation for x(t+h). Since this approximation
may contain x(t + h) on both sides, we call this method an implicit method. In contrast, Euler's
Method, and the Runge-Kutta methods in the following section are explicit methods: they can be
used to directly compute an approximate step. We make this clear by examples.
Example 10.4. Consider the ODE
dx(t)
~St = cosx >
x(0) = 2vr.
In this case, if we wish to approximate x(0 + h) using Backwards Euler's Method, we have to find
x{h) such that x(h) = x(0) + cosx(h). This is equivalent to finding a root of the equation
g(y) = 2vr + cosy - y =
The function g(y) is nonincreasing (g'(y) is nonpositive), and has a zero between 6 and 8, as seen
in Figure 10.4. The techniques of Chapter 4 are needed to find the zero of this function.
Figure 10.4: The nonincreasing function g(y) = 2tt + cosy — y is shown.
10.2. RUNGE-KUTTA METHODS
143
Example 10.5. Consider the ODE from Example 10.1:
dx(t)
dt
x(0)
■■ X,
1.
The actual solution is x(t) = e t . Euler's Method underestimates x(t), as shown in Figure 10.1.
Using Backwards Euler's Method, gives the approximation:
x(t + h) = x(t) + hx(t + h),
x(t)
x(t + h)
1 - h
Using Backwards Euler's Method gives an overestimate, as shown in Figure 10.5. This is because
each step proceeds on the tangent line to a curve ke t to a point on that curve. Generally Backwards
Euler's Method is preferred over vanilla Euler's Method because it gives equivalent stability for
larger stepsize h. This effect is not evident in this case.
Figure 10.5: Backwards Euler's Method applied to approximate x' = x, x(0) = 1. The approxi-
mation is an overestimate, as each step is to a point on a curve ke t , through the tangent at that
point. Compare this figure to Figure 10.1, which shows the same problem approximated by Euler's
Method, with the same stepsize, h = 0.3.
10.2 Runge-Kutta Methods
Recall the ODE problem: find some x(t) such that
f(t,x(t)),
dx(t)
dt
144
CHAPTER 10. ORDINARY DIFFERENTIAL EQUATIONS
where /, a, c are given.
Recall that the Taylor's series method has problems: either you are using the first order method
(Euler's Method), which suffers inaccuracy, or you are using a higher order method which requires
evaluations of higher order derivatives of x(t), i.e., x"(t), x'"(t), etc. This makes these methods less
useful for the general setting of / being a blackbox function. The Runge-Kutta Methods (don't
ask me how it's pronounced) seek to resolve this.
10.2.1 Taylor's Series Redux
We fall back on Taylor's series, in this case the 2-dimensional version:
°° 1 / 8 8\ i
f(* + h,y + k) = ^-(h- + k-) f(x,y).
The thing in the middle is an operator on f(x,y). The partial derivative operators are interpreted
exactly as if they were algebraic terms, that is:
K'Tx +k Ty) I(X ' V)
There is a truncated version of Taylors series:
/(, + h,y + k) - g g [k Tx + k-j /(«,,) 4- (» s + k-j /(«,«.
where (x, y) is a point on the line segment with endpoints (x, y) and (x + h, y + k).
For practice, we will write this for n = 1:
f(x + h,y + k) = f(x, y) + hf x (x, y) + kf y (x, y) + - (h 2 f xx (x, y) + 2hkf xy (x, y) + k 2 f yy (x, y))
10.2.2 Deriving the Runge-Kutta Methods
We now introduce the Runge-Kutta Method for stepping from x(t) to x(t + h). We suppose that
a, f3,w\,W2 are fixed constants. We then compute:
K X i-hf(t,x)
K 2 <- hf(t + ah,x + pKi}.
Then we approximate:
x(i + h)<r- x(t) + wi^Ti + w 2 K 2 .
We now examine the "proper" choices of the constants. Note that w\ = l,w 2 = corresponds
to Euler's Method, and does not require computation of K 2 . We will pick another choice.
= f(x,y),
h df(x,y) h df(x,y)
dx
dy
h 2d ^m^ + 2 hk ^im + k 2d2f ^
dx 2
dxdy
dy 2
10.2. RUNGE-KUTTA METHODS
145
Also notice that the definition of K 2 should be related to Taylor's theorem in two dimensions.
Let's look at it:
K 2 /h = f(t + ah,x + PK 1 )
= f(t + ah,x + (3hf(t,x))
= f + ahf t + /3hff x + \ (a 2 h 2 f tt (i, x) + af3fh 2 f tx (i, x) + 2 h 2 f 2 f x x{i, x)) .
Now reconsider our step:
x(t + h)= x(t) + Wl K x + w 2 K 2
= x(t) + Wl hf + w 2 hf + w 2 ah 2 f t + w 2 (3h 2 ff x + O (h 3 ) .
= x(t) + ( Wl + w 2 ) hx'(t) + w 2 h 2 (aft + pff x ) + O (h 3 ) .
= x(t) + ( Wl + w 2 ) hx'(t) + w 2 h 2 { a f t + Pff x ) + O (h 3 ) .
If we happen to choose our constants such that
wi + w 2 = 1, aw 2 = - = 0w 2 ,
then we get
x(t + h) = x(t) + hx'(t) + -h 2 {f t + ff x ) + O (h 3 )
= x(t) + hx'(t) + hi 2 x"{t) + O (h 3 ) ,
i.e., our choice of the constants makes the approximate x(t + h) good up to a O (/i 3 ) term, because
we end up with the Taylor's series expansion up to that term. cool.
The usual choice of constants is a = (5 = 1, w\ = w 2 = \. This gives the second order Runge-
Kutta Method :
x(t + h)^ x(t) + -f{t, x) + -/(/ + h, x + hf(t, x)).
This can be written (and evaluated) as
K!<-hf(t,x)
K 2 *- hf(t + h,x + Ki)
x(t + h)<- x(t) + - (Kx + K 2 )
Another choice is a = (3 = 2/3, w\ = l/4,w 2 = 3/4. This gives
h 3/i
x (t + h)^x(t) + -f(t,x) + —f[l
(10.4)
The Runge-Kutta Method of order two has error term O {h 3 ). Sometimes this is not enough
and higher-order Runge-Kutta Methods are used. The next Runge-Kutta Method is the order four
method:
<-hf(t,x)
K 2
«- hf{t + h/2,x + K 1 /2)
K 3
^hf(t + h/2,x + K 2 /2)
K 4
^hf(t + h,x + K 3 )
x(t + h)
x(t) + \(Ki + 2K 2 + 2K 3 + K 4 ) .
(10.5)
146
CHAPTER 10. ORDINARY DIFFERENTIAL EQUATIONS
This method has order O (/i 5 ) . (See http : //mathworld. wolf ram. com/Runge-KuttaMethod.html)
The Runge-Kutta Method can be extrapolated to even higher orders. However, the number
of function evaluations grows faster than the accuracy of the method. Thus the methods of order
higher than four are normally not used.
We already saw that w\ = 1, W2 = corresponds to Euler's Method. We consider now the case
wi = 0, W2 = 1. The method becomes
X (t + h )^ x ( t ) + hf(t+^,x+ lf(t, s)) .
This is called the modified Euler's Method . Note this gives a different value than if Euler's
Method was applied twice with step size h/2.
10.2.3 Examples
Example 10.6. Consider the ODE:
f x' = (te) s -(f) 2
\ x(l) = 1
Use h = 0.1 to compute x(l.l) using both Taylor's Series Methods and Runge-Kutta methods
of order 2.
10.3 Systems of ODEs
Recall the regular ODE problem: find some x(t) such that
^ = f(t,x(t)),
x(a) = c,
where /, a, c are given.
Sometimes the physical systems we are considering are more complex. For example, we might
be interested in the system of ODEs:
^l = f(t,x(t),y(t)),
*M=g(t,x(t),y(t)),
x(a) = c,
y(a) = d.
Example 10.7. Consider the following system of ODEs:
f dx(t)
dt
t 2 - X
dt — y t y l,
x(0) = 1,
{ y(o) = o.
You should immediately notice that this is not a system at all, but rather a collection of two ODEs:
dx(t)
dt
x(0) = 1,
t 2 - X
10.3. SYSTEMS OF ODES
147
and
Ml -v 2 + v
dt — y ^ y
2/(0) = 0.
These two ODEs can be solved separately. We call such a system uncoupled. In an uncoupled
system the function f(t,x,y) is independent of y, and g(t,x,y) is independent of x. A system
which is not uncoupled, is, of course, coupled. We will not consider uncoupled systems.
10.3.1 Larger Systems
There is no need to stop at two functions. We may imagine we have to solve the following problem:
find x\ (t) , X2 (t) , . . . , x n {t) such that
f ^ = fi(t, Xl (t),x 2 (t),...,x n (t)),
<^ = f 2 (t,x 1 (t),x 2 (t),...,x n (t)),
^■ = fn{t,X 1 {t),X 2 (t),...,X n (t)),
xi(a) = ci, x 2 (a) = c 2 , . . .,x n (a) = c n .
The idea is to not be afraid of the notation, write everything as vectors, then do exactly the
same thing as for the one dimensional case! That's right, there is nothing new but notation:
Let
Xl
- x' } -
' fi '
Cl
X =
X2
, x' =
x 2
, F =
h
, c =
C2
. X n .
fn
Cn
Then we want to find X such that
X'(t) = F(t,X(t))
X (a) = C.
(10.6)
Compare this to equation 10.1.
10.3.2 Recasting Single ODE Methods
We will consider stepping methods for solving higher order ODEs. That is, we know X(t), and
want to find X (t + h) . Most of the methods we looked at for solving one dimensional ODEs can
be rewritten to solve higher dimensional problems.
For example, Euler's Method, equation 10.3 can be written as
X(t + h)^ X(t) + hF (t, X(t)) .
As in ID, this method is derived from approximating X by its linearization.
In fact, our general strategy of using Taylor's Theorem also carries through without change.
That is we can use the k th order method to step as follows:
X (t + h)^ x(t) + hx'(t) + |x"(t) + . . . + ^x {h) (t).
148 CHAPTER 10. ORDINARY DIFFERENTIAL EQUATIONS
Additionally we can write the Runge-Kutta Methods as follows:
Order two:
<- hF (t,X)
K-2
^hF(t + h,X + K 1 )
X(t + h)
^X(t) + ±(K 1 + K 2 ).
Order four:
K x <- hF(t,X)
K 2 <— hF (^ t + h,X + ^K 1
K s <— hF ( t +±h,X + ^K 2
K 4 ^hF(t + h,X + K 3 )
X(t + h)^ x(t) + i (JTi + 2K 2 + 2K 3 + K 4 ) .
Example Problem 10.8. Consider the system of ODEs:
f x[(t) =xi-x\
x' 2 (t) = t + X\ + x 3
x' 3 (t) = x 2 — x\
*i(0) = 1,
x 2 (0) = 0,
I x 3 (0) = 1.
Approximate X(0.1) by taking a single step of size 0.1, for Euler's Method, and the Runge-Kutta
Method of order 2.
Solution: We write
x 2 (t) - xl(t)
" 1 "
F{t,X(t)) =
t + Xi(t) + x 3 {t)
X(0) =
x 2 {t)-x\{t)
1
Thus we have
F(0,X(0))
Euler's Method then makes the approximation
X(0.1) <- X(0) + O.LF(0, X(0))
The Runge-Kutta Method computes:
K x <- 0.1F(0,X(0)) =
-1
2
-1
" 1 "
" -0.1 "
" 0.9 "
+
0.2
0.2
1
-0.1
0.9
-0.1
0.2
-0.1
and K 2 ^0.1F(0A,X(0) + K 1 )
-0.061
0.19
-0.061
10.3. SYSTEMS OF ODES
149
Then
X(0.1)^X(0) + -(K 1 + K 2 )
' 1 "
" -0.0805 "
" 0.9195 "
+
0.195
0.195
1
-0.0805
0.9195
10.3.3 It's Only Systems
The fact we can simply solve systems of ODEs using old methods is good news. It allows us to
solve higher order ODEs. For example, consider the following problem: given f,a,co,ci, . . . , c n ,
find x(t) such that
^(t) = /(^(t),x'(f),..,xM)(i)),
x(a) = cq, x'(a) = ci, x"(a) = c 2 , . . . , x^ n ~ 1 \a) =
c n -i-
We can put this in terms of a system by letting
X q — X ^ 3^ ]_ — «^ 5 X 2 — ^ • • • ^ — 1 — *^ ^ ^ •
Then the ODE plus these n — 1 equations can be written as
a4_i(*) = /(*,»()(*), zi(i),x 2 (<),..., x„_i(i))
*{,(*)= Xi(t)
xi(t) = i 2 (<)
<- 2 (*) = X n -l(t)
, x (a) = c , xi(a) = ci, x 2 (a) = c 2 , . . . , x n _i(a) = c n _i.
This is just a system of ODEs that we can solve like any other.
Note that this trick also works on systems of higher order ODEs. That is, we can transform
any system of higher order ODEs into a (larger) system of first order ODEs.
Example Problem 10.9. Rewrite the following system as a system of first order ODEs:
x'>(t) = (t/y(t))+x>(t)-l
y (*) = (x'(t)+»(t))
I x(0) = 1, x'(0) = 0,2/(0) = -1
Solution: We let xq = x, x± = x', x 2 = y, then we can rewrite as
f a^t) = (t/x 2 (t))+x 1 (t) - 1
4(t) = xi(t)
X 2(*) = ( Xl {t)+x 2 (t))
{ x (Q) = 1, xi(0) = 0,x 2 (0) =
150
CHAPTER 10. ORDINARY DIFFERENTIAL EQUATIONS
10.3.4 It's Only Autonomous Systems
In fact, we can use a similar trick to get rid of the time-dependence of an ODE. Consider the first
order system
fi (t,X 1 ,X 2 , ...,x n ),
f 2 (t,x 1 ,x 2 ,...,x n ),
f 3 (t,Xl,X 2 , ...,x n ),
X n — f n (t,X±, X 2 , . . . , X n ) ,
xi(a) = ci, x 2 (a) = c 2 , . . .,x n (a) = c n .
We can get rid of the time dependence and thus make the system autonomous by making t
another variable function to be found. We do this by letting xq = t, then adding the equations
x'q = 1, and xo(a) = a, to get the system:
x'l = fl (X ,X!,X 2 ,
x' 2 = f 2 (x ,xi,x 2 ,
x' 3 = h {X ,X!,X 2 ,
X n — f n (xq, X\, X 2 , . . . , Xn) ,
k x (a) = a, xi(a) = ci, x 2 (a) = c 2 , . . .,x n (a) = c n .
An autonomous system can be written in the more elegant form:
X' = F(X)
X (a) = C.
Moreover, we can consider the phase curve of an autonomous system.
One can consider X to be the position of a particle which moves through R n over time. For
an autonomous system of ODEs, the movement of the particle is dependent only on its current
position and not on time. 1 A phase curve is a flow line in the vector field F. You can think of
the vector field F as being the velocity, at a given point, of a river, the flow of which is time
independent. Under this analogy, the phase curve is the path of a vessel placed in the river. In our
deterministic world, phase curves never cross. This is because at the point of crossing, there would
have to be two different values of the tangent of the curve, i.e., the vector field would have to have
two different values at that point.
The following example illustrates the idea of phase curves for autonomous ODE.
Example 10.10. Consider the autonomous system of ODEs, without an initial value:
x[(t) = -2(x 2 -2.4)
x' 2 (t) = 3(xi - 1.2)
We can rewrite this ODE as
X'(t) = F(X(t))
1 Note that if you have converted a system of n equations with a time dependence to an autonomous system, then
the autonomous system should be considered a particle moving through R n+1 . In this case time becomes one of the
dimensions, and thus we speak of position, instead of time.
10.3. SYSTEMS OF ODES
This is an autonomous ODE. The associated vector field can be expressed as
151
F(X) =
This vector field is plotted in Figure 10.6.
-2
3
4.8
-3.6
Figure 10.6: The vector field F (X) is shown in R 2 . The tips of the vectors are shown as circles.
(Due to budget restrictions, arrowheads were not available for this figure.)
You should verify that this ODE is solved by
X(t) =
V2cos(VQt + t )
" 1.2 "
+
2.4
for any r > 0, and to £ (0, 2ir] . Given an initial value, the parameters r and to are uniquely
determined. Thus the trajectory of X over time is that of an ellipse in the plane. 2 Thus the
family of such ellipses, taking all r > 0, form the phase curves of this ODE. We would expect the
approximation of an ODE to never cross a phase curve. This is because the phase curve represents
the
Euler's Method and the second order Runge-Kutta Method were used to approximate the ODE
with initial value
" 1.5
X(0)
1.5
2 In the case r = 0, the ellipse is the "trivial" ellipse which consists of a point.
152
CHAPTER 10. ORDINARY DIFFERENTIAL EQUATIONS
Figure 10.7: The simulated phase curves of the system of Example 10.10 are shown for (a) Euler's
Method and (b) the Runge-Kutta Method of order 2. The Runge-Kutta Method used larger step
size, but was more accurate. The actual solution to the ODE is an ellipse, thus the "spiralling"
observed for Euler's Method is spurious.
10.3. SYSTEMS OF ODES
153
Euler's Method was used for step size h = 0.005 for 3000 steps. The Runge-Kutta Method was
employed with step size h = 0.015 for 3000 steps. The approximations are shown in Figure 10.7.
Given that the actual solution of this initial value problem is an ellipse, we see that Euler's
Method performed rather poorly, spiralling out from the ellipse. This must be the case for any step
size, since Euler's Method steps in the direction tangent to the actual solution; thus every step of
Euler's Method puts the approximation on an ellipse of larger radius. Smaller stepsize minimizes
this effect, but at greater computational cost.
The Runge-Kutta Method performs much better, and for larger step size. At the given resolu-
tion, no spiralling is observable. Thus the Runge-Kutta Method outperforms Euler's Method, and
at lesser computational cost. 3
3 Because the Runge-Kutta Method of order two requires two evaluations of F , whereas Euler's Method requires
only one, per step, it is expected that they would be 'evenly matched' when Runge-Kutta Method uses a step size
twice that of Euler's Method.
154
CHAPTER 10. ORDINARY DIFFERENTIAL EQUATIONS
Exercises
(10.1) Given the ODE:
r x '(t) = t 2 - x 2
\ x(l) = 1
Approximate x(l.l) by using one step of Euler's Method. Approximate x(l.l) using one
step of the second order Runge-Kutta Method.
(10.2) Given the ODE:
x[(t) = X\X 2 + t
x' 2 (t) = 2X2 — x 2
xi(l) =
I x 2 (l) = 3
Approximate X(l.l) by using one step of Euler's Method. Approximate (1.1) using one
step of the second order Runge-Kutta Method.
(10.3) Consider the separable ODE
x'(t) = x(t)g(t),
for some function g. Use backward's Euler (or, alternatively, the right endpoint rule for
approximating integrals) to derive the approximation scheme
/ , \ x (t)
J < t + fc) ^ i-MtW
What could go wrong with using such a scheme?
(10.4) Consider the ODE
x '{t) = x{t)+t 2 .
Using Taylor's Theorem, derive a higher order approximation scheme to find x(t + h) based
on x(t), t, and h.
(10.5) Which of the following ODEs can you show are stable? Which are unstable? Which seem
ambiguous?
(a) x'(t) = —x — arctanx (b) x'(t) = 2x + tan a; (c) x'(t) = —Ax — e x (d) x'(t) = x + x 3
(e) x'it) = t + cos x
(10.6) Rewrite the following higher order ODE as a system of first order ODEs:
' x"(t) = x(t)- sin £'(*),
< x(0) = 0,
k X'(0) = 7T.
(10.7) Rewrite the system of ODEs as a system of first order ODEs
' x"(t)=y(t)+t-x(t),
y'(t) = x'{t) + x(t) - 4,
< x'(0) = 1,
x(0) = 2,
I y(0) = 0.
(10.8) Rewrite the following higher order system of ODEs as a system of first order ODEs:
f x"'(t)=y"(t)-x'(t),
y"'(t)=x"(t)+y'(t),
X (0) = x"(0) = y'(0) = -1,
k x'(0) = y(0) = y"(0) = 1.
10.3. SYSTEMS OF ODES
155
(10.9) Implement Euler's Method for approximating the solution to
= /(*,*(*)),
x(a) = c.
Your m-file should have header line like:
function xfin = euler (f , a, c ,h, steps)
where xfin should be the approximate solution to the ODE at time a + steps *h.
(a) Run your code with f (t, x) = — x, using a = / c, and using h * steps = 1.0. In this
case the actual solution is
xfin = ce" 1 .
This ODE is stable, so your approximation should be good (cf. Example 10.3). Exper-
iment with different values of h.
(b) Run your code with f(t,x) = x, using a = / c, and using h * steps = 1.0. In this
case the actual solution is
xfin = ce.
Your actual solution may be a rather poor approximation. Prepare a plot of error
versus h for varying values of h. (cf. Figure 7.1 and Figure 10.2)
(c) Run your code with f (t, x) = 1/ (1 — x) , using a = 0, c = 0.1, and using h*steps = 2.0.
Plot your approximations for varying values of steps. For example, try steps values of
35, 36, and 90, 91. Also try large values of steps, like 100, 500, 1000. Can you guess what
the actual solution is supposed to look like? Can you explain the poor performance of
Euler's Method for this ODE?
(10.10) Implement the Runge-Kutta Method of order two for approximating the solution to
= /(*,*(*)),
x(a) = c.
Your m-file should have header line like:
function xfin = rkm(f , a, c ,h, steps)
where xfin should be the approximate solution to the ODE at time a + steps *h.
Test your code with the ODEs from the previous problem. Does the Runge-Kutta Method
perform any better than Euler's Method for the last ODE?
CHAPTER 10. ORDINARY DIFFERENTIAL EQUATIONS
Appendix A
Old Exams
A.l First Midterm, Fall 2003
PI (7 pnts) Clearly state Taylor's Theorem for f(x + h). Include all hypotheses, and state the
conditions on the variable that appears in the error term.
P2 (7 pnts) Write out the Lagrange Polynomial £o{x) for nodes xo,xi, . . . ,x n . That is, write the
polynomial that has value 1 at xq and has value at x\, X2, ■ ■ ■ , x n .
P3 (7 pnts) Write the Lagrange form of the polynomial of lowest degree that interpolates the
following values:
Xk
-1
1
Vk
11
17
P4 (14 pnts) Use the method of Bisection for finding the zero of the function f(x) = x 3 — 3x 2 +
5x — |. Let ao = 0, and &o = 2. What are 02, 62?
P5 (21 pnts) • Write out the iteration step of Newton's Method for finding a zero of the func-
tion f(x). Your iteration should look like
x k +i = ??
• Write the Newton's Method iteration for finding y/b.
• Letting xq = 1, find the x\ and X2 for your iteration. Note that y/5 ~ 2.236.
P6 (7 pnts) Write out the iteration step of the Secant Method for finding a zero of the function
f(x). Your iteration should look like
x k+ i =??
P7 (21 pnts) Consider the following approximation:
, ^ /(x + fe)-5/(x) + 4/(x+|)
3h
• Derive this approximation using Taylor's Theorem.
• Assuming that f{x) has bounded derivatives, give the accuracy of the above approxima-
tion. Your answer should be something like O
• Let f(x) = x 2 . Approximate /'(0) with this approximation, using h = |.
P8 (21 pnts) Consider the data:
k
1
2
Xk
3
1
-1
f(x k )
5
15
9
157
158
APPENDIX A. OLD EXAMS
• Construct the divided differences table for the data.
• Construct the Newton form of the polynomial p(x) of lowest degree that interpolates
f(x) at these nodes.
• Suppose that these data were generated by the function f(x) = x 3 — 5x 2 + 2x + 17.
Find a numerical upper bound on the error \p(x) — f{x)\ over the interval [—1,3]. Your
answer should be a number.
A. 2 Second Midterm, Fall 2003
PI (6 pnts) Write the Trapezoid rule approximation for f(x)dx, using the nodes {xi}™ =0 .
P2 (6 pnts) Suppose Gaussian Elimination with scaled partial pivoting is applied to the matrix
0.0001
-0.1
-0.01
0.001
43
91
-43
1
-12
26
-1
-7
l
2
-10
7T
Which row contains the first pivot? You must show your work.
P3 (14 pnts) Let cf)(n) be some quadrature rule that divides an interval into 2 n equal subintervals.
Suppose
f f(x) dx = <j>(n) + a 5 h 5 n + a w h™ + a 15 h£ + ...,
J a
where h n = Using evaluations of (/>(•), devise a "Romb erg-style" quadrature rule that is
a O approximation to the integral.
P4 (14 pnts) Compute the LU factorization of the matrix
A =
1 -1
2
2 2-1
using naive Gaussian Elimination. Show all your work. Your final answer should be two
matrices, L and U. Verify your answer, i.e., verify that A = LU.
P5 (14 pnts) Run one iteration of either Jacobi or Gauss Seidel on the system:
1
3
4 "
' -3 "
-2
2
1
x =
-3
-1
3
-2
1
Start with = [1 1 1] T . Clearly indicate your answer x^ > . Show your work. You must
indicate which of the two methods you are using.
P6 (21 pnts) Derive the two-node Chebyshev Quadrature rule:
J 1 J{x)dx^c[f(x {i ) + f(x 1 )}.
Make this rule exact for all polynomials of degree < 2.
P7 (28 pnts) Consider the quadrature rule:
f 1 f(x) dx « zof{0) + zi/(l) + z 2 f'(0) + z 3 f'(l).
Jo
A.3. FINAL EXAM, FALL 2003
159
• Write four linear equations involving Zi to make this rule exact for polynomials of the
highest possible degree.
• Solve your system for z using naive Gaussian Elimination.
A. 3 Final Exam, Fall 2003
PI (4 pnts) Clearly state Taylor's Theorem for f(x + h). Include all hypotheses, and state the
conditions on the variable that appears in the error term.
P2 (4 pnts) State the conditions that characterize S(x) as a natural cubic spline on the interval
[a, b]. Let the knots be a = to < t\ < ti < . . . < t n = b.
P3 (4 pnts) Let xq = 3,x± = — 5,X2 = 0. Write out the Lagrange Polynomial ^o(^) fo r these
nodes; that is, write the polynomial that has value 1 at xo and has value at xi,X2-
P4 (4 pnts) Consider the ODE:
x'(t) = f(t,x(t))
x{a) = c
Write out the Euler Method to step from x(t) to x(t + h).
P5 (7 pnts) Compute the LU factorization of the matrix
A =
-1 1
1 1 -1
-2 2
using naive Gaussian Elimination. Show all your work. Your final answer should be two
matrices, L and U. Verify your answer, i.e., verify that A = LU.
P6 (9 pnts) Let a be given. Determine a quadrature rule of the form
J —t
f{x) dx « A/(0) + Bf"(-a) + Cf"(a)
that is exact for polynomials of highest possible degree. What is the highest degree polynomial
for which this rule is exact?
P7 (9 pnts) Suppose that 4>{n) is some computational function that approximates a desired quan-
tity, L. Furthermore suppose that
4>(n)
_1 _2 _3
L + a\n 2 + a2fi 2 + a^n 2 + . . .
Combine two evaluations of </>(■) in the manner of Richardson to get aO(n x ) approximation
to L.
P8 (9 pnts) Find the least-squares best approximate of the form y = In (cx) to the data
x
Xl
2/o
V\
Vn
Caution: this should not be done with basis vectors.
P9 (9 pnts) Let
T = {/(x) = CoX + Cl logX + C2 COSX I Co, C\,C2 6l).
Set up (but do not attempt to solve) the normal equations to find the least-squares best / G T
to approximate the data:
x
Xl
2/o
2/1
2/n
160
APPENDIX A. OLD EXAMS
The normal equations should involve the vector c = [co c\ C2] , which uniquely determines a
function in T .
P10 (9 pnts) Consider the following approximation:
/"(*)
3/(x - fe) - 4/(x) + /(x + 3h)
6h 2
• Derive this approximation using Taylor's Theorem.
• Assuming that /(x) has bounded derivatives, give the accuracy of the above approxima-
tion. Your answer should be something like O
Pll (11 pnts) Consider the ODE:
x'(t) = x{t)g(t)
x(0) = 1
Use the Fundamental Theorem of Calculus to write a step update of the form:
rt+h
:{t + h) = x(t) + J ??
• Approximate the integral using the Trapezoid Rule to get an explicit step update of the
form
x(t + h) «- Q(t,x{t),h)
• Suppose that g(t) > m > for all t. Do you see any problem with this stepping method
iih> 2/m?
Hint: the actual solution to the ODE is
P12 (9 pnts) Consider the ODE:
x(t) = e& 9 ^ Ar .
x'{t) = e x ^ +sm(x(t))
x(0) =
Write out the Taylor's series method of order three to approximate x(t + h) given x(t). It is
permissible to write the update as a small program like:
At) ^ Qo(x(t))
X"(t) <r- Q^xit^At))
x(t + h) «- R(x(t),x'{t),x"{t),...)
rather than write the update as a monstrous equation of the form
x{t + h) *- M(x(t)).
P13 (15 pnts) Consider the data:
k
1
2
Xk
-0.5
0.5
f(Xk)
5
15
9
• Construct the divided differences table for the data.
A.4. FIRST MIDTERM, FALL 2004
161
• Construct the Newton form of the polynomial p(x) of lowest degree that interpolates
f(x) at these nodes.
• Suppose that these data were generated by the function f(x) = 40x 3 — 32x 2 — 6x + 15.
Find a numerical upper bound on the error \p(x) — f(x)\ over the interval [—0.5,0.5].
Your answer should be a number.
P14 (15 pnts) • Write out the iteration step of Newton's Method for finding a zero of the
function f(x). Your iteration should look like
Xk+l =??
• Write the Newton's Method iteration for finding ^/25/4.
• Letting xq = |, find the x\ and X2 for your iteration. Note that ^25/4 ~ 1.842.
P15 (15 pnts) • Let x be some given number, and let /(•) be a given, "black box" function.
Using Taylor's Theorem, derive a O (h) approximation to the derivative f'{x). Call your
approximation <fi(h); the approximation should be defined in terms of /(•) evaluated at
some values, and should probably involve x and h.
• Let f(x) = x 2 + 1. Approximate /'(0) using your c/)(h), for h=\.
• For the same f(x) and h, get a O (/i 2 ) approximation to /'(0) using the method of
Richardson's Extrapolation; your answer should probably include a table like this:
£>(0,0)
£(1,0) £(1,1)
A.4 First Midterm, Fall 2004
[Definitions] Answer no more than two of the following.
PI (12 pnts) Clearly state Taylor's Theorem for f(x + h). Include all hypotheses, and state the
conditions on the variable that appears in the error term.
P2 (12 pnts) Write the general form of the iterated solution to the problem Ax = b. Let Q be
your "splitting matrix," and use the factor u. Your solution should look something like:
??a; (fe) =??a; (fc - 1) +?? or x^ =??a; (fe - 1) +??
For one of the following variants, describe what Q is, in relation to A : Richardson's, Jacobi,
Gauss-Seidel.
P3 (12 pnts) Write the iteration for Newton's Method or the Secant Method for finding a root to
the equation f(x) = 0. Your solution should look something like:
??
Xk+l =?? + ^
[Problems] Answer no more than four of the following.
PI (14 pnts) Perform bisection on the function graphed below. Let the first interval be [0,256].
Give the second, third, fourth, and fifth intervals.
P2 (14 pnts) Give a O (/i 4 ) approximation to sin (h + 7r/2) . Find a reasonable C such that the
truncation error is no greater than Ch 4 .
P3 (14 pnts) Use Newton's Method to find a good approximation to \/24. Use xq = 3. That is,
iteratively define a sequence with xq = 3, such that Xk — ►
P4 (14 pnts) One of the roots of x 2 — bx + c = as found by the quadratic formula is subject
to subtractive cancellation when \b\ 3> |c| . Which root is it? Rewrite the expression for that
root to eliminate subtractive cancellation.
P5 (14 pnts) Find the LU decomposition of the following matrix by Gaussian Elimination with
naive pivoting.
"3 2 4
-6 1 -6
12 -12 10 _
P6 (14 pnts) Consider the linear equation:
" 6
1
1 "
5 "
2
4
x =
-6
_ 1
2
6
3
Letting x^ = [11 1] T , find a^ 1 ) using either the Gauss-Seidel or Jacobi iterative schemes.
Use weighting u = 1. You must indicate which scheme you are using. Show all your work.
[Questions] Answer only one of the following.
PI (20 pnts) Suppose that the eigenvalues of A are in [a,(3] with < a < (3. Find the to that
minimizes ]|l — u;A 2 || 2 . What is the minimal value of ||l — <^A 2 || 2 for this optimal uil For full
credit you need to make some argument why this cj minimizes this norm.
P2 (20 pnts) Let f(x) = have a unique root r. Recall that the form of the error for Newton's
Method is
p _ -/"(&) ,
efc+1 - m^) k '
where e& = r — x^. Suppose that f(x) has the properties that f(x) > m > for all x, and
\f"( x )\ ^ M for all x. Find 5 > such that je^j < 5 insures that x n — ► r. Justify your answer.
P3 (20 pnts) Suppose that f(r) = = f'(r) / f"(r), and f"(x) is continuous. Letting = x^ — r,
show that for Newton's Method,
1
when |efc| is sufficiently small. (Hint: Use Taylor's Theorem twice.) Is this faster or slower
convergence than for Newton's Method with a simple root?
A. 5 Second Midterm, Fall 2004
[Definitions] Answer no more than two of the following three.
PI (10 pnts) For a set of distinct nodes xo,xi, . . . ,xi3, what is the Lagrange Polynomial ^(a;)?
You may write your answer as a product. What degree is this polynomial?
A.5. SECOND MIDTERM, FALL 2004
163
P2 (10 pnts) Complete this statement of the polynomial interpolation error theorem: Let p be the
polynomial of degree at most n interpolating function f at the n+1 distinct nodes Xo,x±, . . . ,x n
on [a,b]. Let /( ra+1 ) be continuous. Then for each x G [a, b] there is some £ G [a, b] such that
P3 (10 pnts) State the conditions which define the function S(x) as a natural cubic spline on the
interval [a, b] .
[Problems] Answer no more than five of the following six.
PI (16 pnts) How large must n be to interpolate the function e x to within 0.001 on the interval
[0, 2] with the polynomial interpolant at n equally spaced nodes?
P2 (16 pnts) Let 4>(h) be some computable approximation to the quantity L such that
,15
4>(h) = L + a 5 h b + ai /i + a i5 h
Combine evaluations of </>(•) to devise a O (h w ) approximation to L.
P3 (16 pnts) Derive the constants A and B which make the quadrature rule
j:
f(x) dx « Af
15
+ Bf(Q) + Af
T5
exact for polynomials of degree < 2. Use this quadrature rule to approximate
2
f 1
TT = —
J-l 1
+ X 2
■ dx
P4 (16 pnts) Find the polynomial of degree < 3 interpolating the data
X
-l
1
2
y
-2
8
4
(Hint: using Lagrange Polynomials will probably take too long.)
P5 (16 pnts) Find constants a, (3 such that the following is a quadratic spline on [0,4].
Q(x) = <
x 2 - 3x + 4
a(x - l) 2 + f3(x - 1) + 2
,2
^x z + x — 4
P6 (16 pnts) Consider the following approximation:
f( x + h)-5f(x) + 4f(x+l)
< x < 1
1 < x < 3
3 < x < 4
3h
• Derive this approximation via Taylor's Theorem.
• Assuming that f(x) has bounded derivatives, give the accuracy of the above approxima-
tion. Your answer should be something like O (h 7 ).
• Let f(x) = x 2 . Approximate f'(0) with this approximation, using h = |.
164
APPENDIX A. OLD EXAMS
A.6 Final Exam, Fall 2004
[Definitions] Answer all of the following.
PI (10 pnts) Clearly state Taylor's Theorem for f(x + h). Include all hypotheses, and state the
conditions on the variable that appears in the error term.
P2 (10 pnts) Write the iteration for Newton's Method or the Secant Method for finding a root to
the equation f(x) = 0. Your solution should look something like:
Xk+l =?? +
??
77
P3 (10 pnts) Write the general form of the iterated solution to the problem Ax = b. Let Q be
your "splitting matrix," and use the factor ui. Your solution should look something like:
77 x ( k ) =??a;( A;_1 )+??
or
x (k) = ?? iC ( fc - 1 )_)_??
For one of the following variants, describe what Q is, in relation to A : Richardson's, Jacobi,
Gauss-Seidel.
P4 (10 pnts) Define what it means for the function f(x) € T to be the least squares best approx-
imant to the data
Xq
Xi
yo
yi
Vn
[Problems] Answer all of the following.
PI (10 pnts) Rewrite this system of higher order ODEs as a system of first order ODEs:
x"(t)
= x'{t) [x(t) + ty{t)\
y"(t)
= y'(t) \y{t) - tx(t)}
s'(0)
= 1
< x(0) =
= 4
1/(0)
= -3
y(o) --
= -2
P2 (10 pnts) Consider the ODE:
V(t) =t(x + l) 2
x(2) = 1/2
Use Euler's Method to find an approximate value of x(2.1) using a single step.
P3 (10 pnts) Let (j>{h) be some computable approximation to the quantity L such that
4>{h) = L + a 4 /i 4 + a 8 /i 8 + a 12 h lz + ...
Combine evaluations of <p(-) to devise a O (/i 8 ) approximation to L.
P4 (15 pnts) Consider the approximation
/'(*)
9 f{x + h)- 8f(x) - f{x - 3h)
12h
(a) Derive this approximation using Taylor's Theorem.
(b) Assuming that f{x) has bounded derivatives, give the accuracy of the above approxima-
tion. Your answer should be something like O (h 7 ).
A.6. FINAL EXAM, FALL 2004
165
P5 (15 pnts) Find constants, a, (3, 7 such that the following is a natural cubic spline on [1,3].
Q(x)
'a (x - l) 3 + (x - l) 2 + 12 (x - 1) + 1 : 1 < x < 2
7 x d + 18x^ - 30x + 19
: 2 < x < 3
P6 (15 pnts) Devise a subroutine that, given an input number z, calculates an approximate to
via Newton's Method. Write your subroutine in pseudo code, Matlab, or C/C++. Include
reasonable termination conditions so your subroutine will terminate if it finds a good approx-
imate solution. Your subroutine should use only the basic arithmetic operations of addition,
subtraction, multiplication, division; in particular your subroutine should not have access to
a cubed root or logarithm operator.
P7 (15 pnts) Consider the ODE:
jx'(t) = cos (x(t)) + sin (x(t))
\x(0) = 1
Use Taylor's theorem to devise a O (/i 3 ) approximant to x(t + h) which is a function of x, t,
and h only. You may write your answer as a program:
x >(t) <- Q (x(t))
a/'(t)<-Q 1 (x(t),x'(t))
x <- R(x{t),x'{t),x"(t),...)
such that x(t + h) = x + O (h 3 ) .
[Questions] Answer all of the following.
PI (20 pnts) Consider the "quadrature" rule:
/' f(x) dx « A f(0) + Ai/(1) + A 2 f'(l) + A 3 f"(l)
Jo
(a) Write four linear equations involving A^ to make this rule exact for polynomials of the
highest possible degree.
(b) Find the constants Ai using Gaussian Elimination with naive pivoting.
P2 (20 pnts) Consider the data:
k
1
2
x k
0.25
0.5
f(x k )
1.5
1
0.5
(a) Construct the divided differences table for the data.
(b) Construct the Newton form of the polynomial p(x) of lowest degree that interpolates
/(x) at these nodes.
(c) What degree is your polynomial? Is this strange?
(d) Suppose that these data were generated by the function
f(x) = 1 +
COS 2"7TX
Find an upper bound on the error \p(x) — f(x)\ over the interval [0,0.5]. Your answer
should be a number.
166
APPENDIX A. OLD EXAMS
P3 (20 pnts) Let
T = {f(x) = coVx + ci£(x) + c 2 e x | c , c±, c 2 eK},
where £(x) is some black box function. We are concerned with finding the least-squares best
/ G T to approximate the data
x
Xl
Xn
yo
yi
Vn
(a) Set up (but do not attempt to solve) the normal equations to find the the vector c =
[co ci C2] T , which determines the least squares approximant in T .
(b) Now suppose that the function l(x) happens to be the Lagrange Polynomial associated
with X7, among the values XQ,X\,...,x n . That is, £(xj) = 5jj. Prove that f(xj) = j/7,
where / is the least squares best approximant to the data.
P4 (10 pnts) State some substantive question which you thought might appear on this exam, but
did not. Answer this question (correctly).
Appendix B
GNU Free Documentation License
Version 1.2, November 2002
Copyright ©2000,2001,2002 Free Software Foundation, Inc.
59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Everyone is permitted to copy and distribute verbatim copies of this license document, but changing it is not allowed.
Preamble
The purpose of this License is to make a manual, textbook, or other functional and useful document "free" in the sense of
freedom: to assure everyone the effective freedom to copy and redistribute it, with or without modifying it, either commercially
or noncommercially. Secondarily, this License preserves for the author and publisher a way to get credit for their work, while
not being considered responsible for modifications made by others.
This License is a kind of " copy left" , which means that derivative works of the document must themselves be free in the
same sense. It complements the GNU General Public License, which is a copyleft license designed for free software.
We have designed this License in order to use it for manuals for free software, because free software needs free documentation:
a free program should come with manuals providing the same freedoms that the software does. But this License is not limited
to software manuals; it can be used for any textual work, regardless of subject matter or whether it is published as a printed
book. We recommend this License principally for works whose purpose is instruction or reference.
1. APPLICABILITY AND DEFINITIONS
This License applies to any manual or other work, in any medium, that contains a notice placed by the copyright holder
saying it can be distributed under the terms of this License. Such a notice grants a world-wide, royalty-free license, unlimited
in duration, to use that work under the conditions stated herein. The "Document", below, refers to any such manual or work.
Any member of the public is a licensee, and is addressed as "you". You accept the license if you copy, modify or distribute the
work in a way requiring permission under copyright law.
A "Modified Version" of the Document means any work containing the Document or a portion of it, cither copied
verbatim, or with modifications and/or translated into another language.
A "Secondary Section" is a named appendix or a front-matter section of the Document that deals exclusively with the
relationship of the publishers or authors of the Document to the Document's overall subject (or to related matters) and contains
nothing that could fall directly within that overall subject. (Thus, if the Document is in part a textbook of mathematics, a
Secondary Section may not explain any mathematics.) The relationship could be a matter of historical connection with the
subject or with related matters, or of legal, commercial, philosophical, ethical or political position regarding them.
The "Invariant Sections" are certain Secondary Sections whose titles are designated, as being those of Invariant Sections,
in the notice that says that the Document is released under this License. If a section does not fit the above definition of Secondary
then it is not allowed to be designated as Invariant. The Document may contain zero Invariant Sections. If the Document does
not identify any Invariant Sections then there are none.
The "Cover Texts" are certain short passages of text that are listed, as Front-Cover Texts or Back-Cover Texts, in
the notice that says that the Document is released under this License. A Front-Cover Text may be at most 5 words, and a
Back-Cover Text may be at most 25 words.
A "Transparent" copy of the Document means a machine-readable copy, represented in a format whose specification is
available to the general public, that is suitable for revising the document straightforwardly with generic text editors or (for
images composed of pixels) generic paint programs or (for drawings) some widely available drawing editor, and that is suitable
for input to text formatters or for automatic translation to a variety of formats suitable for input to text formatters. A copy
made in an otherwise Transparent file format whose markup, or absence of markup, has been arranged to thwart or discourage
subsequent modification by readers is not Transparent. An image format is not Transparent if used for any substantial amount
of text. A copy that is not "Transparent" is called "Opaque".
167
168
APPENDIX B. GNU FREE DOCUMENTATION LICENSE
Examples of suitable formats for Transparent copies include plain ASCII without markup, Texinfo input format, LaTeX
input format, SGML or XML using a publicly available DTD, and standard-conforming simple HTML, PostScript or PDF
designed for human modification. Examples of transparent image formats include PNG, XCF and JPG. Opaque formats
include proprietary formats that can be read and edited only by proprietary word processors, SGML or XML for which the
DTD and/or processing tools are not generally available, and the machine-generated HTML, PostScript or PDF produced by
some word processors for output purposes only.
The "Title Page" means, for a printed book, the title page itself, plus such following pages as are needed to hold, legibly,
the material this License requires to appear in the title page. For works in formats which do not have any title page as such,
"Title Page" means the text near the most prominent appearance of the work's title, preceding the beginning of the body of
the text.
A section "Entitled XYZ" means a named subunit of the Document whose title either is precisely XYZ or contains XYZ
in parentheses following text that translates XYZ in another language. (Here XYZ stands for a specific section name mentioned
below, such as "Acknowledgements", "Dedications", "Endorsements", or "History".) To "Preserve the Title" of
such a section when you modify the Document means that it remains a section "Entitled XYZ" according to this definition.
The Document may include Warranty Disclaimers next to the notice which states that this License applies to the Document.
These Warranty Disclaimers are considered to be included by reference in this License, but only as regards disclaiming warranties:
any other implication that these Warranty Disclaimers may have is void and has no effect on the meaning of this License.
2. VERBATIM COPYING
You may copy and distribute the Document in any medium, cither commercially or noncommercially, provided that this
License, the copyright notices, and the license notice saying this License applies to the Document are reproduced in all copies,
and that you add no other conditions whatsoever to those of this License. You may not use technical measures to obstruct or
control the reading or further copying of the copies you make or distribute. However, you may accept compensation in exchange
for copies. If you distribute a large enough number of copies you must also follow the conditions in section 3.
You may also lend copies, under the same conditions stated above, and you may publicly display copies.
3. COPYING IN QUANTITY
If you publish printed copies (or copies in media that commonly have printed covers) of the Document, numbering more
than 100, and the Document's license notice requires Cover Texts, you must enclose the copies in covers that carry, clearly and
legibly, all these Cover Texts: Front-Cover Texts on the front cover, and Back-Cover Texts on the back cover. Both covers must
also clearly and legibly identify you as the publisher of these copies. The front cover must present the full title with all words of
the title equally prominent and visible. You may add other material on the covers in addition. Copying with changes limited to
the covers, as long as they preserve the title of the Document and satisfy these conditions, can be treated as verbatim copying
in other respects.
If the required texts for either cover are too voluminous to fit legibly, you should put the first ones listed (as many as fit
reasonably) on the actual cover, and continue the rest onto adjacent pages.
If you publish or distribute Opaque copies of the Document numbering more than 100, you must either include a machine-
readable Transparent copy along with each Opaque copy, or state in or with each Opaque copy a computer-network location from
which the general network-using public has access to download using public-standard network protocols a complete Transparent
copy of the Document, free of added material. If you use the latter option, you must take reasonably prudent steps, when you
begin distribution of Opaque copies in quantity, to ensure that this Transparent copy will remain thus accessible at the stated
location until at least one year after the last time you distribute an Opaque copy (directly or through your agents or retailers)
of that edition to the public.
It is requested, but not required, that you contact the authors of the Document well before redistributing any large number
of copies, to give them a chance to provide you with an updated version of the Document.
4. MODIFICATIONS
You may copy and distribute a Modified Version of the Document under the conditions of sections 2 and 3 above, provided
that you release the Modified Version under precisely this License, with the Modified Version filling the role of the Document,
thus licensing distribution and modification of the Modified Version to whoever possesses a copy of it. In addition, you must
do these things in the Modified Version:
A. Use in the Title Page (and on the covers, if any) a title distinct from that of the Document, and from those of previous
versions (which should, if there were any, be listed in the History section of the Document). You may use the same title
as a previous version if the original publisher of that version gives permission.
B. List on the Title Page, as authors, one or more persons or entities responsible for authorship of the modifications in the
Modified Version, together with at least five of the principal authors of the Document (all of its principal authors, if it
has fewer than five), unless they release you from this requirement.
C. State on the Title page the name of the publisher of the Modified Version, as the publisher.
D. Preserve all the copyright notices of the Document.
E. Add an appropriate copyright notice for your modifications adjacent to the other copyright notices.
169
F. Include, immediately after the copyright notices, a license notice giving the public permission to use the Modified Version
under the terms of this License, in the form shown in the Addendum below.
G. Preserve in that license notice the full lists of Invariant Sections and required Cover Texts given in the Document's
license notice.
H. Include an unaltered copy of this License.
I. Preserve the section Entitled "History", Preserve its Title, and add to it an item stating at least the title, year, new
authors, and publisher of the Modified Version as given on the Title Page. If there is no section Entitled "History" in
the Document, create one stating the title, year, authors, and publisher of the Document as given on its Title Page,
then add an item describing the Modified Version as stated in the previous sentence.
J. Preserve the network location, if any, given in the Document for public access to a Transparent copy of the Document,
and likewise the network locations given in the Document for previous versions it was based on. These may be placed
in the "History" section. You may omit a network location for a work that was published at least four years before the
Document itself, or if the original publisher of the version it refers to gives permission.
K. For any section Entitled "Acknowledgements" or "Dedications", Preserve the Title of the section, and preserve in the
section all the substance and tone of each of the contributor acknowledgements and/or dedications given therein.
L. Preserve all the Invariant Sections of the Document, unaltered in their text and in their titles. Section numbers or the
equivalent are not considered part of the section titles.
M. Delete any section Entitled "Endorsements". Such a section may not be included in the Modified Version.
N. Do not retitle any existing section to be Entitled "Endorsements" or to conflict in title with any Invariant Section.
O. Preserve any Warranty Disclaimers.
If the Modified Version includes new front-matter sections or appendices that qualify as Secondary Sections and contain
no material copied from the Document, you may at your option designate some or all of these sections as invariant. To do this,
add their titles to the list of Invariant Sections in the Modified Version's license notice. These titles must be distinct from any
other section titles.
You may add a section Entitled " Endorsements" , provided it contains nothing but endorsements of your Modified Version by
various parties-for example, statements of peer review or that the text has been approved by an organization as the authoritative
definition of a standard.
You may add a passage of up to five words as a Front-Cover Text, and a passage of up to 25 words as a Back-Cover Text,
to the end of the list of Cover Texts in the Modified Version. Only one passage of Front-Cover Text and one of Back-Cover
Text may be added by (or through arrangements made by) any one entity. If the Document already includes a cover text for
the same cover, previously added by you or by arrangement made by the same entity you are acting on behalf of, you may not
add another; but you may replace the old one, on explicit permission from the previous publisher that added the old one.
The author(s) and publisher(s) of the Document do not by this License give permission to use their names for publicity
for or to assert or imply endorsement of any Modified Version.
5. COMBINING DOCUMENTS
You may combine the Document with other documents released under this License, under the terms defined in section 4
above for modified versions, provided that you include in the combination all of the Invariant Sections of all of the original
documents, unmodified, and list them all as Invariant Sections of your combined work in its license notice, and that you preserve
all their Warranty Disclaimers.
The combined work need only contain one copy of this License, and multiple identical Invariant Sections may be replaced
with a single copy. If there arc multiple Invariant Sections with the same name but different contents, make the title of each
such section unique by adding at the end of it, in parentheses, the name of the original author or publisher of that section if
known, or else a unique number. Make the same adjustment to the section titles in the list of Invariant Sections in the license
notice of the combined work.
In the combination, you must combine any sections Entitled "History" in the various original documents, forming one section
Entitled " History" ; likewise combine any sections Entitled " Acknowledgements" , and any sections Entitled " Dedications" . You
must delete all sections Entitled "Endorsements".
6. COLLECTIONS OF DOCUMENTS
You may make a collection consisting of the Document and other documents released under this License, and replace the
individual copies of this License in the various documents with a single copy that is included in the collection, provided that
you follow the rules of this License for verbatim copying of each of the documents in all other respects.
You may extract a single document from such a collection, and distribute it individually under this License, provided you
insert a copy of this License into the extracted document, and follow this License in all other respects regarding verbatim
copying of that document.
7. AGGREGATION WITH INDEPENDENT WORKS
170
APPENDIX B. GNU FREE DOCUMENTATION LICENSE
A compilation of the Document or its derivatives with other separate and independent documents or works, in or on a
volume of a storage or distribution medium, is called an "aggregate" if the copyright resulting from the compilation is not used
to limit the legal rights of the compilation's users beyond what the individual works permit. When the Document is included
in an aggregate, this License does not apply to the other works in the aggregate which are not themselves derivative works of
the Document.
If the Cover Text requirement of section 3 is applicable to these copies of the Document, then if the Document is less than
one half of the entire aggregate, the Document's Cover Texts may be placed on covers that bracket the Document within the
aggregate, or the electronic equivalent of covers if the Document is in electronic form. Otherwise they must appear on printed
covers that bracket the whole aggregate.
8. TRANSLATION
Translation is considered a kind of modification, so you may distribute translations of the Document under the terms of
section 4. Replacing Invariant Sections with translations requires special permission from their copyright holders, but you may
include translations of some or all Invariant Sections in addition to the original versions of these Invariant Sections. You may
include a translation of this License, and all the license notices in the Document, and any Warranty Disclaimers, provided that
you also include the original English version of this License and the original versions of those notices and disclaimers. In case
of a disagreement between the translation and the original version of this License or a notice or disclaimer, the original version
will prevail.
If a section in the Document is Entitled "Acknowledgements", "Dedications", or "History", the requirement (section 4) to
Preserve its Title (section 1) will typically require changing the actual title.
9. TERMINATION
You may not copy, modify, sublicense, or distribute the Document except as expressly provided for under this License.
Any other attempt to copy, modify, sublicense or distribute the Document is void, and will automatically terminate your rights
under this License. However, parties who have received copies, or rights, from you under this License will not have their licenses
terminated so long as such parties remain in full compliance.
10. FUTURE REVISIONS OF THIS LICENSE
The Free Software Foundation may publish new, revised versions of the GNU Free Documentation License from time to
time. Such new versions will be similar in spirit to the present version, but may differ in detail to address new problems or
concerns. See http://www.gnu.org/copylcft/.
Each version of the License is given a distinguishing version number. If the Document specifies that a particular numbered
version of this License "or any later version" applies to it, you have the option of following the terms and conditions either of
that specified version or of any later version that has been published (not as a draft) by the Free Software Foundation. If the
Document does not specify a version number of this License, you may choose any version ever published (not as a draft) by the
Free Software Foundation.
ADDENDUM: How to use this License for your documents
To use this License in a document you have written, include a copy of the License in the document and put the following
copyright and license notices just after the title page:
Copyright ©YEAR YOUR NAME. Permission is granted to copy, distribute and/or modify this document
under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free
Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. A copy of
the license is included in the section entitled " GNU Free Documentation License" .
If you have Invariant Sections, Front-Cover Texts and Back-Cover Texts, replace the "with. ..Texts." line with this:
with the Invariant Sections being LIST THEIR TITLES, with the Front-Cover Texts being LIST, and with the
Back-Cover Texts being LIST.
If you have Invariant Sections without Cover Texts, or some other combination of the three, merge those two alternatives
to suit the situation.
If your document contains nontrivial examples of program code, we recommend releasing these examples in parallel under
your choice of free software license, such as the GNU General Public License, to permit their use in free software.
Bibliography
[1] Guillaume Bal. Lecture notes, course on numerical analysis, 2002. URL
http : //www . Columbia . edu/~gb2030/C0URSES/E6302/NumAnal . ps.
[2] Samuel Daniel Conte and Carl de Boor. Elementary Numerical Analysis; An Algorithmic Ap-
proach. McGraw-Hill, New York, third edition, 1980. ISBN 0070124477.
[3] James F. Epperson. An Introduction to Numerical Methods and Analysis. Wiley, 2001. ISBN
0-471-31647-4.
[4] Richard Hamming. Numerical Methods for Scientists and Engineers. Dover, New York, 1987.
ISBN 0486652416.
[5] Eugene Isaacson and Herbert Bishop Keller. Analysis of Numerical Methods. Dover, New York,
1966. ISBN 0486680290.
[6] David Kincaid and Ward Cheney. Numerical analysis. Brooks/Cole Publishing Co., Pacific
Grove, CA, second edition, 1996. ISBN 0-534-33892-5. Mathematics of scientific computing.
[7] David Kincaid and Ward Cheney. Numerical Mathematics and Computing. Brooks/Cole Pub-
lishing Co., Pacific Grove, CA, fourth edition, 1999. ISBN 0-534-35184-0.
[8] James Stewart. Calculus, Early Trans cendentals. Brooks/Cole-Thomson Learning, Belmont,
CA, fifth edition, 2003. ISBN 0-534-39330-6.
171
Index
Backwards Euler's Method, 141
big-O, 3
bisection method, 50
centered difference, 91
Chebyshev Polynomials, 69, 122
Chebyshev Quadrature, 114
composite quadrature rule, 100
Dirichlet function, 99
divided differences, 67
eigenvalue, 8
eigenvector, 8
elementary row operations, 25
error
roundoff, 89, 138
truncation, 89, 138
Euler's Method, 136, 137
Gaussian Elimination, Naive, 25
Gaussian Quadrature, 107
Heaviside function, 98
Hilbert Matrix, 48
inner product, 6
Intermediate Value Theorem, 49
knot, 79
Lagrange Polynomials, 63
least squares, 119
definition, 117
linear, 119
LU Factorization, 33
method of undetermined coefficients, 108
multiplier, 28
natural cubic spline, 83
Newton's Method, 52
norm, 7
172
subordinate, 9
normal equations, 118
ODE
errors, 138
higher order methods, 137
stability, 138
stepping, 136
partition, 79
phase curve, 150
Richardson Extrapolation, 92, 104
Riemann Integrable, 98
Riemann integral, 97
Romberg Algorithm, 104, 105
Runge Function, 69
Runge-Kutta Method, 143
secant method, 58
splines, 79
basis splines, 84
linear, 79
quadratic, 81
subtractive cancellation, 4
Taylor's Series Method, 136
Taylor's Theorem, 1, 3
trapezoidal rule, 100
error, 102
recursive, 106
vector space, 5
Live Music
Archive
Librivox
Free Audio
Metropolitan Museum
Cleveland
Museum of Art
Internet
Arcade
Console Living Room
Open Library
American
Libraries
TV News
Understanding
9/11